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Abstract. In this paper higher order mimetic discretizations are introduced which are firmly 
rooted in the geometry in which the variables are defined. The paper shows how basic constructs 
in differential geometry have a discrete counterpart in algebraic topology. Generic maps which 
switch between the continuous differential forms and discrete cochains will be discussed and 
^ finally a realization of these ideas in terms of mimetic spectral elements is presented, based on 

^ projections for which operations at the finite dimensional level commute with operations at the 

continuous level. 

The two types of orientation (inner- and outer-orientation) will be introduced at the continu- 
QQ ous level, the discrete level and the preservation of orientation will be demonstrated for the new 

I mimetic operators. The one-to-one correspondence between the continuous formulation and the 

discrete algebraic topological setting, provides a characterization of the oriented discrete bound- 
I I ary of the domain. The Hodge decomposition at the continuous, discrete and finite dimensional 

level will be presented. It appears to be a main ingredient of the structure in this framework. 
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1. Introduction 



1.1. Motivation. The starting point of physics is the science of measurable, quantifiable objects 
and the relation among these objects. Any measurable quantity is associated with a spatial and a 
temporal geometric object. For instance, the measurement of velocity in air flow can be performed 
by particle image velocimetry (PIV), where tracer particles are released in the flow, two pictures 
are taken at consecutive time instants and the average velocity is calculated from the distance a 
particle has traveled divided by the time interval between the two snapshots. Velocity, measured 
in this way, is therefore associated with a curve in space (the trajectory of the particle between 
the two time instants) and a time interval (the time interval between the two snapshots). This 
association of velocity with a curve and a time interval does not depend on the particular way in 
which the measurement of velocity is performed. Note that 



is an exact relation between the velocity in the time interval and the position of the tracer 

particles at and i^, irrespective of the particular path traced out by the particle. So, if we could 
measure these two positions with infinite accuracy, we would have the time integral of the velocity 
correct. The approximation enters into the measurement by assuming that the particle moves at 
a constant velocity in a straight line from r(t^) to r(t^), which allows us to equate the velocity in 
the interval [t^ , t^] to 



We can decompose this velocity measurement in two consecutive steps: 1.) A reduction step, 
where we sample the position of the tracer particles at discrete time instants 2.) A reconstruction 
map, where we assume a certain behaviour of the particles between the two time instants. 

These two steps, reduction and reconstruction, which implicitly play an important role in any 
measurement, will also turn out to be the two key ingredients in setting up mimetic discretizations. 

Strictly speaking, the velocity measured in the PIV experiment is only the time-averaged veloc- 
ity, but by assuming that the trajectory of the tracer particle is sufficiently smooth as a function 
of time, we can reduce the time interval such that we can quite accurately determine 'the velocity' 
at a given position and at a given time instant. An alternative approach would be to sample at 
a larger number of time instants, say t^, . . . ,1"" and reconstruct the trajectory based on the mea- 
sured positions. We will call the first approach (reducing the time interval) h-refinement whereas 
the second approach (reconstruction of the trajectory using more time instants) will be called 
p-refinement. Similar ideas are used in discretization where refinement of the mesh is denoted by 
/i- refinement while a reconstruction based on more samples is referred to as p-refinement. 

Ultimately, in many physical theories, one takes the limit for all lengths and time intervals to 
zero, which enables physicists and engineers to talk about the velocity in a point at a certain time 
instant, 'v(t,x,y, z). Any connection with a distance and a time interval is lost after this limiting 
process. Another well-known example is mass contained in a volume, V. The average density is 
the mass divided by the volume. By taking the limit for V — > we obtain the density in a point, 
p{t,x,y, z). Again, the connection with the volume is lost after taking this limit. 

Bear in mind that this limiting process is purely mathematical. If we consider the PIV experi- 
ment again to measure the local velocity, we always need a finite time interval in order to evaluate 
the average velocity. If we would reduce the time interval to zero, no velocity measurement could 
be made. So despite the fact that we can accurately determine the velocity in a flow at a certain 
location and at a certain time, this measured velocity will always be connected to a time interval 
and the displacement along a curve. 

This association of physical variables with spatial and temporal geometric elements can be 
done for all physical variables. This is, however, beyond the scope of this paper, but the interested 
reader is referred to the work of [HI HZl [75] and especially the forthcoming book by Tonti, [50] . 
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Once we acknowledge that there is such an association between physical variables and geometric 
objects, wc need to take orientation of the geometric object into account. A curve with cndpoints 
A and B possesses two orientations. Either the curve from A ^ B is taken as the positively 
oriented curve, or the curve B ^ A is taken to be positively oriented. The same holds for the 
orientation of surfaces (oriented clockwise or counter-clockwise) and three dimensional volumes 
(right-hand rule or left-hand rule). Note that the notion of orientation does not prompt itself 
when we only consider physical variables defined at time instants and in points, although it is 
useful to consider the orientation of points as well. 

Physical variables are associated to geometric objects and geometric objects have an orientation, 
however, the physical quantity is independent of the orientation of the associated geometric object. 
If we choose the direction of time to be positive when pointing in the past, the integral value of 
the velocity changes sign 




So integral values (the ones that are measurable) arc intimately connected to the orientation of 
geometry, while average values - and in the limit densities - are insensitive to the orientation of 
space and time, because 

— / vd,t = — / V dt . 

There also exist physical variables which do change when the orientation is reversed. Therefore, 
we need to distinguish between two types of orientation: inner- and outer orientation. With 
inner orientation, we mean the orientation in the geometric object such as for instance the electric 
current in a wire or the rotation in a plane, whereas outer refers to the orientation outside the 
geometric object such as the Biot-Savart law around the wire and or the flux through the plane. 

It turns out that the association with oriented geometric objects is a vital ingredient in the 
description of physics and when we perform the limiting process to define all physical quantities 
in points and at time instants without reference to the associated geometric objects much of this 
rich structure of the physical model will be lost. 

In this paper, therefore, we want to set up a framework in which we mimic the association 
of physical variables with oriented geometric objects for computational analysis. The aim is to 
develop families of numerical discretizations which work on general quadrilateral grids of arbitrary 
order. By 'families' we mean that we formulate the basic requirements a numerical discretization 
needs to possess in order to be compatible with its associated geometry. This leads, among others, 
to exact discrete representations for the gradient operator, the curl operator and the divergence 
operator. Also, the explicit distinction between inner- and outer-orientation with their associated 
cell complexes, will anatomy of the boundary of the domain. This, in turn, will clarify the issue 
of where and how to prescribe boundary values. 

The mimetic structure that will be introduced in this paper ensures that the reduced physical 
model behaves in the same way as the full infinite dimensional system. Let A and B be two 
physical quantities and T a continuous operator which maps A onto B, T{A) = B, then the 
following diagram commutes 



A - 


B 


TT 






A 


h - 


— ^ Bh 



Here is tt a suitably constructed projection operator which maps continuous variables in finite 
dimensional representations. So tt o T = T o tt, i.e. we can perform the operation T at the 
continuous level and then discretize or first discretize and then apply the operator T. In that 
sense, operations at the discrete level truly mimic the behaviour of the operators at the continuous 
level. 

These properties have quite some consequences for practical applications, but we have chosen 
not to present an extensive gallery of applications. Only a few simple examples will be given which 



MIMETIC FRAMEWORK 



5 



serve to illustrate some of the claims in this paper. Applications of ideas presented in this paper 
can be found in [TUl 113 BH El IMl [3g . 

1.2. Prior and related work. Over the years numerical analysts have developed numerical 
schemes which preserve some of the structure of the differential models they aim to approximate, 
so in that respect the whole mimetic idea is not new. A recent development is that the proper 
language in which to encode these structures/symmetries is the language of differential geometry. 
Another novel aspect of mimetic discretizations is the identification of the metric-free part of 
differential models, which can be conveniently described in terms of algebraic topology which 
employs the strong analogy between differential geometry and algebraic topology. 

The relation between differential geometry and algebraic topology in physical theories was 
first established by Tonti, [75] • Around the same time Dodziuk, [5D], set up a finite difference 
framework for harmonic functions based on Hodge theory. Both Tonti and Dodziuk introduce 
differential forms and cochain spaces as the building blocks for their theory. The relation between 
differential forms and cochains is established by the Whitney map (/c-cochains — > fc-forms) and 
the De Rham map (fc-forms -> fc-cochains). The interpolation of cochains to differential forms 
on a triangular grid was already established by Whitney, [82] . These interpolatory forms are now 
known as the Whitney forms. 

Hyman and Scovel, i33| . set up the discrete framework in terms of cochains, which are the 
natural building blocks of finite volume methods. Later Bochev and Hyman, |6j extended this 
work and derived discrete operators such as the discrete wedge product, the discrete codifferential, 
the discrete inner product, etc. These operators are all cochain operators. 

In a finite difference/ volume context Robidoux, Hyman, Steinberg and Shashkov, [5^H551 [551 [Ml 
1731 1771 178| used symmetry considerations to discretize diffusion problems on rough grids and with 
non-smooth non-isotropic diffusion coefficients. In a recent paper by Robidoux and Steinberg |71j a 
discrete vector calculus in a finite difference setting is presented. It satisfies the discrete differential 
operators grad, curl and div exactly and the numerical approximations are all contained in the 
constitutive relations, which are already polluted by modeling and experimental error. This paper 
also contains an extensive list of references to mimetic methods. For mimetic finite differences, 
see also Brezzi et al. , [TTJ IH] . 

The application of mimetic ideas to unstructured staggered grids has been extensively studied 
by Perot, [58l l60H62| [83] . Especially the recent paper, [59 , lucidly describes the rationale of 
preserving symmetries in numerical algorithms. 

Mattiussi, puts the geometric ideas proposed by Tonti in an finite volume, finite dif- 

ference and finite element context. The idea of switching between cochains and differential forms 
is also prominent in the work of Hiptmair, for instance |30j . This work also displays the close 
connection between finite volume methods and finite element methods. 

Mimetic methods show a clear connection between the variables (differential forms) and the 
geometry in which these variables are defined. The most 'geometric approach' is described in the 
work by Desbrun et al., [T51 IHH |S7| and the thesis by Hirani, IJT]. 

The 'Japanese papers' by Bossavit, [7|[Sjj serve as an excellent introduction and motivation 
for the use of differential forms in the description of physics and the use in numerical modeling. 
The field of application is electromagnetism, but these papers are sufficiently general to extend all 
concepts to other fields of expertise. 

In a series of papers by Arnold, Falk and Winther, [U |31|S], a finite element exterior calculus 
framework is developed. Just like in this paper, Arnold, Falk and Winther consider methods 
of arbitrary order. Higher order methods are also described by Rapetti, ,^03 HI] and Hiptmair, 
[29] . Possible extensions to spectral methods were described by Robidoux, [70] and applications 
of mimetic spectral methods can be found in [TUl HZl US HH |M1 [55] . 

Isogeometric reconstruction was used by Buffa et al., [13], Evans, [H] and Hiemstra, [551 164j . 

Although the cited literature is far from complete, the above references serve as excellent in- 
troduction into the field of mimetic discretization techniques. 

1.3. Scope and outline of this paper. In the introduction and work cited above it has been 
revealed that geometry plays an important role in mimetic methods. In computational engineering 
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one usually works with fields and densities, i.e. the variables obtained after the limiting process. 
The main reason fields have emerged as the preferred way of encoding physics is because physical 
laws can then be stated in terms of differential equations. 

An alternative description is in terms of integral equations. The appeal of an integral approach 
lies in the fact that the physical laws can be stated without any limiting process involved, ren- 
dering them closer to the physical measurement process and more suited for a discrete treatment. 
Integration can be interpreted as duality pairing between geometry and variables connected to 
this geometry (differential forms). 

An important differential operator in differential geometry is the exterior derivative. The exte- 
rior derivative can be defined in terms of geometric concepts, i.e. the boundary operator, through 
the generalized Stokes Theorem. When written in terms of conventional vector calculus, the ex- 
terior derivative is either the gradient, the curl or the divergence, depending on the context. The 
introduction of the exterior derivative allows one to uniquely decompose the space of differen- 
tial forms into a direct sum of sub-spaces. This Hodge decomposition generalizes the classical 
Helmholtz decomposition for non-contractible domains. So, if we want to incorporate geometry 
into our physical description, differential geometry is a concise and potent way to do so. There- 
fore, in Section [2] a brief introduction into differential geometry will be presented. Although this 
material can be found in any book on differential geometry, [T21 [Ml US] , we include this section to 
introduce our notation and in the remainder of the paper we want to highlight which properties 
from differential geometry are retained at the discrete level. 

Despite the fact that physics requires metric concepts like length, angles and area, many struc- 
tures in physics are completely independent of the metric. These non-metric concepts are called 
topological. In Section [3] elements from algebraic topology will be discussed which is required for 
the development of the mimetic spectral element framework. It will be shown that the struc- 
ture of algebraic topology resembles the structure of differential geometry and therefore algebraic 
topology could serve as the discrete setting for our numerical framework. 

In Section |4] the connection between algebraic topology and differential geometry will be estab- 
lished. Based on the existence of a suitable reduction map, TZ, which maps fc-forms onto fc-cochains 
and a reconstruction map, I, which converts /c-cochains to /c-forms, a general mimetic framework 
will be set up using the projection operator tt/i = IoTZ^ which maps the space of differential forms, 
A*^, to a finite dimensional space of differential forms, A^. This section resembles the paper by 
Bochev and Hyman, f6|, but the main difference is that in Section [4] the finite dimensional discrete 
space consists of differential forms while Bochev and Hyman take the cochains as their discrete 
variables. 

In Section [5] the actual polynomial reduction and reconstruction maps are presented which sat- 
isfy the requirements described in Section [4j Their composition forms a bounded linear projection 
as proven in this section. This section is accompanied with many examples of the actions of the 
various operators. 

In Section [6] we will review the tools developed in this paper and look back to the introduction 
and see how this approach may enable us to faithfully simulate problems in physical sciences. 
Furthermore, potential future directions will be identified. 

Although this outline suggests a collection of seemingly unrelated scientific fields, several 
ideas/concepts permeate throughout the paper. Ultimately, all concepts contribute to mimetic, 
numerical concepts: 



(1) The exterior derivative — > coboundary operator — > discrete gradient, curl and diver- 
gence; 

(2) The Hodge decomposition — > cohomology group — > discrete Helmholtz decomposition; 

(3) The wedge product — > tensor products — >■ basis functions on quadrilateral elements; 

(4) The behaviour under mappings: the puUback operator — > the cochain map — > mimetic 
discretization on highly deformed meshes; 

(5) Inner- and outer orientation — > the double De Rham complex — > boundary of the domain 
— > the trace operator/boundary values — > spectral method on a staggered grid; 



MIMETIC FRAMEWORK 



7 



(6) The generalized Stokes Theorem — > discrete generahzed Stokes — > exact conservation 
and existence of scalar and vector potentials. 
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2. Differential Geometric Concepts 



Differential geometry, along with algebraic topology presented in the next section, constitutes 
the basis of the numerical framework presented in this paper. In contrast to conventional vector 
calculus, which is a topic well known to all the readers differential geometry is not a familiar topic 
to most of the readers. Since differential geometry is essential, we include a short introduction 
in order to make this work as much as possible self-contained and to be able to draw analogies 
between differential geometry, algebraic topology and the mimetic scheme we will present. Only 
those concepts from differential geometry which will play a role in the remainder of this article are 
introduced. We start by introducing the concept of manifold, which is the playground in which 
everything is defined, the geometry, and the concept of orientation is presented. Next, differential 
forms, their definition, the operators defined on them and their transformation under mappings 
are introduced subsequently. Differential form spaces will be introduced, including the concept of 
Hodge decomposition. The matter presented here constitutes the mathematical tools with which 
the physical quantities and the physical laws will be represented in the continuous world, and 
what will be mimicked in the numerical framework presented. For a more in depth treatment of 
differential geometry in physics, we refer to 1 ^ [5 ^ [57 1 [501 [721 [75] . 

2.1. Manifolds. The concepts which will be introduced all exist associated to sets endowed with 
enough structure so that one can "do calculus" and which are denoted by manifolds. In M.^ these 
are commonly referred to as points, lines, surfaces and volumes. Generalizing to any dimension a 
manifold can be defined in the following way. 



Definition 1 (Manifold) . '54'1 A /c-dimensional manifold is a set Ai , together with a countable 
collection of subset Ua C Ai, called coordinate charts, and one-to-one functions (pM,a ■ Ua — > Va 
onto connected open subsets Va ofR'^, called local coordinate maps, as in Figure\^ which satisfy 
the following properties: 

(1) The coordinate charts cover Ai: 




Figure 1. Coordinate charts on a manifold. 



a 



(2) On the overlap of any pair of coordinate charts Ua f^Up, the composite map 
is a smooth (infinitely differentiable) function. 
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(3) If X £ lAa and y e Uj^ are distinct points in A4, then there exist open subsets Wa of 
^M.aix) in Va and Wp of ipM./six) in Vp such that 

Since the image of each point p G {Ua H A4) by (pM,a is a point in M*^, it can be written as a k- 
tuple of real numbers: (pM.a{p) = {x^{p), ■ ■ ■ ,x"{p)). This /c-tuple is caUed the local coordinates 
of p and Ua H A4 the coordinate neighborhood. The pair (Ua, (pM.a) is called a local chart or 
local coordinate system. An atlas on a manifold is a collection A = {{Ua, 'i^M,a)} of charts of 
A4, such that IJa^" = -M; the collection of open sets {Ua} constitutes an open covering of the 
manifold M. 

Definition 2 (Maps between manifolds). Let he a k- dimensional smooth manifold and M 
an l-dimensional smooth manifolds. The map $ : — !■ A/" maps between manifolds, if for every 
coordinate chart (pM,a ■ Ua — > Vq C M*^ on Ai and every chart fjs/.p : Up Vp £1 M} on J\f, the 
composite map 

is a smooth map wherever it is defined. See Figure^ for a pictorial representation of the mapping 
between manifolds. 




Figure 2. Mapping between two manifolds, M and Af. 

If the map ^p^f^p o $ o (fij^^ is of maximal rank at p S A4, then there are local coordinates 
X = {x^, . . . , x^) near p and y — (y^ , . . . ,y^) near q = ^{p) G Af such that these coordinates have 
the simple form 

y ^ {x\...,x'',0,...,0) , iil>k, 

or 

y ^ {x\...,x^) , ifl<k. 

Definition 3 (Submanifold). Let M. be a smooth manifold. A submanifold of M. is a subset 
S d M, together with a smooth one-to-one mapping $ : 5 — ^ 5 C satisfying the maximal rank 
condition everywhere, where the parameter space S is some other manifold and S — $(5) is the 
image of In particular, the dimension of S is the same as that of S , and does not exceed the 
dimension of Ai. 
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An important concept is the boundary of a manifold. This concept plays an essential role in 
the generalized Stokes theorem, Theorem [T] to be introduced later on. 

Definition 4 (Complement, interior point, exterior point, boundary point, open set 
and closed set). I37j Given a subset S of a manifold A4, the complement of S in A4 is the set 
of points S'^ := {p € Ai\p ^ S}. Let the ball B^{p) = {x & M \ d{x,p) < e}, then 

(1) A point p is an interior point of S if there exists e > such that the neighborhood (p) 
around p has the property that -Be(p) C S. One writes p G int(5). 

(2) A point p is an exterior point of S if there exists e > such that the neighborhood B^{p) 
around p has the property that Bt{p) n 5 = 0. One writes p G cxt(5). 

(3) A point p is a boundary point of S if every neighborhood B^(p) around p with e > has 
the property that B^{p) H S ^ and B^{p) H 5'^ 7^ 0. One write p € dS . 

A set S is open if, and only if, S = int(5). A set S is closed if, and only if, its complement S"^ 
is open. 

In order to introduce the concept of boundary one needs to introduce the closed upper half-space 
of dimension n: 

H" = {(a;\...,x")eM"|a;">0}, 
with the subspace topology of M" . From Definition |4] it follows that points with > are the 
interior points of H", the points with < are the exterior points of H" and the points with 
x" = are the boundary points of H". 

Proposition 1. \81f Let P and Q be subsets o/H" and ^ : P ^ Q a diffeomorphism. Then $ 
maps interior points to interior points and boundary points to boundary points. 

One can then define an n-manifold with boundary: 

Definition 5 (n-Manifold with boundary, interior point and boundary point of an 
rt- manifold with boundary). ISl^ An n-manifold with boundary, M., is a topological space 
which is locally H". A point p of M is an interior point if there is a chart ^PM,a in which ipM,aip) 
is an interior point o/H". In the same way, a point p is a boundary point of Ai if fM,a{p) is a 
boundary point of IP. The set of boundary points of M is denoted by dAi. 

Definition 6 (Boundary operator). Given an n-manifold with boundary, Ai, the boundary 
operator d is a map d : Ai ^ dAi . 

Corollary 1 (The boundary of a submanifold). Since any submanifold S <Z Ai is a manifold 
in its own right, the boundary of a submanifold is defined as in Definition 

Proposition 2 (Boundary is mapped into a boundary and the boundary is indepen- 
dent of chart). Given two n-dimensional manifolds with boundary, Ai and A" , and a mapping 
(diffeomorphism) between them, $ : — A" , then the interior points and boundary points of Ai 
are mapped onto interior points and boundary points of M, respectively. That is: 

(2.1) d'!>{M) = <^{dM). 

Moreover, interior and boundary points are independent of the choice of chart. 

Proof For the first statement, let tpM.a : Ua d At ^ R'^ and ^pj^/^p : Up d N ^ M.^ , then 
o o fj^ji Q, : M*^ — ^ K*^ is a diffeomorphism and according to Proposition [T] this maps 
boundary points onto boundary points and interior points onto interior points. 

As for the second statement, one takes M = N and in this case the two charts will be '.pM,a 
and fM,i3 and the mapping between the two charts (change of coordinates) will be given by 
Vm,I3 o $ o if'J^ ^ — (pM,i3 ° fMa^ since $ in this case is the identity map, and again ipM,i3 ° a ■ 
M'^ ^ M*^ is a diffeomorphism and according to Proposition [l] this maps boundary points onto 
boundary points and interior points onto interior points. □ 

In n-dimcnsional space it is possible to define n + 1 sub-manifolds of dimension 0, 1, . . . , n, 
respectively. For the case n = 3 one can define, points, lines, surfaces and volumes. Moreover, it 



MIMETIC FRAMEWORK 



11 



is also possible to orient these objects. By orientation one means the generalization of concepts 
such as left and right, front and back, clockwise or counterclockwise, outward and inward, etc. 
The different kinds of orientation will be presented for all the geometric objects that exist in 
3-manifolds (points, lines, surfaces and volumes), as well as generalizations for geometric objects 
of arbitrary dimension. The distinction between inner orientation (solely related to the geometric 
object) and outer orientation (related to both the object and the embedding space) will be given. 
For a more detailed discussion on orientation we recommend [J [121 SZl [751 [7S] . 

The concept of orientation on manifolds is a generalization of the one for vector spaces and 
hence, by extension, to M.^. One starts with the notion of orientation in vector spaces and the 
charts (p will induce an orientation on the manifold 
In 



an orientation is one of the two possible directions, see Figure 3(a) 



In 



is one of the two possible rotations, clockwise or counterclockwise, see Figure 3(b) 



an orientation 
In 



an 



orientation is one of the two possible screw senses, upward clockwise or upward counterclockwise, 



see Figure 3(c) 




(a) ID vector space (b) 2D vector space (c) 3D vector space 

Figure 3. Possible orientations of vector spaces. 

The question is how to generalize this heuristic definition to higher dimensions. In a vector 
space of dimension k one can transform one set of basis vectors, {vi, V2, ■ ■ ■ ,Vk}, into another one, 
{ui,U2, ■ ■ ■ ,Uk}, in the following way: 

Ui = ^ SijVj , 
j 

where Sij are the coefficients of the transformation matrix S with det(S') =/= 0. 

Definition 7 (Orientation). In a vector space of dimension k, orientation is an equivalence 
class of ordered sets of basis vectors whose equivalence relation states that two sets of basis vectors 
belong to the same equivalence class if the transformation matrix, S, between them has det(S') > 0. 

Since the determinant of a change of basis is either positive or negative, there are only two such 
classes. Hence, any vector space has only two orientations. If one chooses one of them, arbitrarily, 
the other one is said to be the opposite orientation. This also corresponds to assigning one of 
them as positive and the other as negative. It is simple to see that this definition of orientation is 
equivalent for the three cases presented above and generalizes this concept to any dimension k. 

As seen in Definition [l] A:- manifolds are locally like TpM*"' = IR'^, and therefore, locally like a 
/c-dimensional vector space. In this way, orientation on a manifold can be defined as: 

Definition 8 (Inner orientation on a manifold). A manifold M is said to be inner oriented 
if for any two overlapping charts, iUa,'-PM.a) and (Up,ipM,p)j of its atlas. A, the Jacobian deter- 
minant of the transformation fj^^ o (/J;vt,Q positive. The inner orientation being the one of the 
equivalence classes of the sets of basis vectors of the tangent space at each point, TpAi, associated 
to these charts. 

If one considers the space in which the fc-manifold A4 is embedded to be M", with dimension 
n > k, then at each point on the manifold there is a space perpendicular to the tangent space, 
TpM, denoted by T^j-Ai whose dimension is {n — k); the normal bundle to the submanifold A4 
embedded in K". Then TpM" — TpM ® T^M. This allows one to define an outer orientation of 
a manifold as: 
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Definition 9 (Outer orientation on a manifold). Consider an oriented k -manifold Ai, with in- 
ner orientation {ui , U2 ,■■ ■ atTpAd, embedded in an n- dimensional Euclidean space. An outer 
orientation of Ai is an orientation for the perpendicular vector space at each point in the manifold, 
T^Ai, {uk+i, Uk+2, ■ ■ ■ , Un} such that all the oriented basis {ui, U2, ■ ■ ■ , Uk, Ufe+i, • ■ ■ , Un} o,re in 
one of the two equivalence classes of the embedding space TpM". 

The particular cases of inner orientation of a O-manifold (point) and outer orientation of a 
n-manifold embedded in a n-dimensional space are treated in a similar way. In both cases, the 
tangent space (points) and the perpendicular space (n-manifolds) have dimension zero. Therefore, 
the points and the 71-manifolds are simply considered as sources or sinks and their orientation can 
be seen as simply induced by the inner orientation of the lines stemming out of them (points) or 
by the outer orientation of it's faces (71-manifold) , see [47l [79] . 

Example 1 (Outer orientation of points). The outer orientation of points depends on the 
embedding space M" . In Figure [J a graphical representation of the outer orientation of a point in 
M" is given for n — 0, ... ,3 




Figure 4. Outer orientation of a point embedded in M", M^, and (from 
left to right). 

Example 2 (Outer orientation of line segments). In Figure^the outer orientation of a line 
segment embedded in M" is shown for n = 1, . . . , 3. 




Figure 5. Outer orientation of a line segment embedded in M^, and (from 
left to right). 

Example 3 (Outer orientation of surfaces). Figure^shows the embedding of a surface in 
M", n = 2 and n = 3. 




Figure 6. Outer orientation of a surface embedded in and M.^ (from left to right). 

Example 4 (Outer orientation of volumes). Figure^shows the outer orientation of a volume 
m R^. 

Figure [S] presents above a sequence of outer-oriented geometric objects of increasing dimension, 
and below a sequence of inner-oriented geometric objects of decreasing order. The objects are 
aranged in such a way that it reveals the similarities with the double de Rham complex and action 
of the Hodge-* operator. Both are introduced in this section later on. 
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Figure 7. Outer orientation of a volume embedded in E'^. 



Outer Orientation 




Figure 8. Inner and outer orientation of /c-manifolds, fc = {0, 1, 2, 3}, in . 



Following the ideas pointed in |80j . it is important to stress that orientation plays an essential 
role in integration. Although not explicitly expressed in the more common differential formulation, 
it appears implicitly by the usage of the right hand rule, for example, or when Stokes' theorem 
is considered, Theorem [T] In this case, the orientation of a fc-manifold and its boundary have to 
be compatible when integrating, for the theorem to hold. Moreover, the distinction between inner 
and outer oriented manifolds will play a central role in the discretization to be presented in this 
paper. This distinction motivates the use of dual grids. 

This work is restricted to manifolds for which a global consistent orientation can be defined, 
called orientable manifolds. 

2.2. Differential forms. Differential forms will play a central role in the development of our 
numerical framework. 

Definition 10 (Differential forms). |2l[53[7^ yl differential k-form, a^^), k > 1 IS a mapping: 

a'-''^ : TpM x • • • x TpM U , 
fe 

which is skew symmetric. That is, for any permutation P of the indices 1, 2, • • • , k: 

a'-'^^vi,--- ,Vk) =sgn(F)a('')(i;p(i),--- ,Vp(^k)) , 

for any ui, • ■ ■ , Wfc € TpAi, see |7y. A 0-form, a^^\ is defined simply as a standard scalar function 
on A4. The space of k-forms on the manifold M is denoted by K^[M). a'-'^-' = when k < or 
k > n. 
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Definition 11 (Wedge product). JJl [75]/ The wedge product, A, of two differential fi 
A''{M) and e A\M) is a mapping A : A'=(7W) x A'(7W) ^ A'=+'(A^), such that: 

(2.2a) (a^'^) + &('') A cS"'^ = a^^^ A c^") + A c*") (Distrihutivity) 

(flC^) A &(')) A c(") = flC^) A (6^') A c^")) 

(2.2b) = a^*^) A A c(™) ('ylssociatiTO^yj 

(2.2c) aaC') A 6^'^ = a^^) A ah^^^ = a(a('=) A 6^')) (Multiplication by scalars) 

(2.2d) a^^) A = (-1)'='6(') A a'^'''' (Skew symmetry) 



orms a 



(k) 



From property (2.2d) we we have 
(2.3) a^'') A a^*^) = 0, Va^'^^ e A^ fc is odd or > ^ . 

Proposition 3. ^ On a manifold Ai of dimension n the space of 1-forms A^(A^) is a linear 
vector space of dimension n that is spanned by n basis elements. A canonical basis, given a local 
coordinate system {x^ , ■ ■ ■ ,a;") is {dx^ , ■ ■ ■ ,dx'^}. The space of k-forms A''(A^) is a linear vector 
space of dimension md has a canonical basis given by: 

{dx'' A • • • A dx'^ll <ii<---<ik<n}. 
Example 5. Examples of 0- forms, 1-forms, 2- forms and 3- forms are given by: 
= a{x, y, z) 

= bi{x,y,z)dx + b2(x,y,z)dy + bz{x,y,z)dz 



(2.4) 



.(2) - 



ci{x, y, z)dyAdz + C2{x, y, z)d2;Ada; + 03(0;, y, z)dxAdy 



w'^^^ — w{x,y, z)dxAdyAdz 

Example 6 (a*^' A a^^^ in M^). dx A dx = dy A dy = dz A dz = . 

Example 7 (a^^) A a^^) in M"*). In R'^, let a^^^ = da; A dy + dz A dt, then 

a(2) A a(2) = (dx A dy + dz A dt) A (dx A dy + dz A di) = 2 dx A dy A dz A dt 7^ . 

Definition 12 (Inner product fc-forms). /7^/ The space of k-forms, A'^(Al), can be equipped 
with a pointwise positive- definite inner product, (•, •) : A^{A4) x A'^(A^) — > M., such that: 



(4i)a...a4^6Wa...a6W) 



"1 ' "fe 



where a\^\b^j^^ £ A^{M), i,j = ,k and (^^i^'? ^j^"*^ '■s inner product of 1-forms, induced 

by the inner product on tangent vectors and by the duality pairing between vectors and 1-forms, 
see m \7&j , and given by: 

where a^^-* — "^^Oidx^, b^^^ = '^jbjdx^ and y*^ are the coefficients of the inverse of the metric 
tensor of rank two, see fJl [75^ . 

2.3. Differential forms under mappings. It is also important to determine how differential 
forms transform under a mapping between two manifolds, see Figure [2] for an example for n = 2 
and where (x, y) = $(^,77), and how integration of fc-forms over manifolds is defined. 

Definition 13 (Fullback operator). JTJ {E^ Consider two manifolds of dimension n, A4 and 
M and a mapping between them, $ : — > Af, such that local coordinates ^* in Ai are mapped 
into local coordinates = $*(^^, • • • ,^") in M . Then the pullback, $*, of a k-form, fc > 1, a^^^ 
is given by: 

a>'^(aW)(t;i,--- ,«fc) :=a(^-)($,(t/i),--- ,<i>.(z;fe)) , 
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where ^i,{v) is the usual pushforward of a vector v, see |7J \24]l - The pullback of a 0-form, a^"-' G 
A°(A/'), is given simply by the composition of the maps 

For a l-form we get the following. 

Proposition 4. If a^^^ £ A^{Af) is given in a local coordinate system as a^^' = ctidx' 

then the pullback $* : K^{M) — ^ A^(A^) is given by: 

(2-5) <J'*(-^^') = E-^Sde^ 

Proposition 5. \32^ The pullback $* has the following properties: 
(2.6a) ^^(ae^) + foC')) = ^^(aC^')) + '^*{b'^''^) (Linearity) 

(2.6b) <^>*{a^^'^ ^b^^^)^<^>*{a^^^)^<^*{b'^^^) (Algebra homomorphism) 

(2.6c) ($2 o •i'l)* = <I'i o ^2 (Composition) 

Integration of differential forms can now be defined in the following way. 



Definition 14 (Pullback integration formulation). \70j The integral of a differential k-form, 
a^^\ over a manifold M. of dimension k is: 

(2.7) / / (^-i)*(aW), 

where (p^^ : S" C M'^ — )• is the inverse of the global chart from A4 to S CZ M.'' . On the right hand 



side of {2.7) one has the usual integral in S d 



Integration can be considered a duality pairing between a differential fc-form and a fc-dimensional 
manifold in the following way 



(2.8) {a^^>,M) / a^ 

JM 

If integration is interpreted as the duality pairing between differential forms and geometry, then 
the relation between a mapping and the associated pullback satisfies 

Proposition 6. Given a mapping $ : Al — > A/", its associated pullback $* and a differential form 
£ A'^(A/') then the following holds: 

(2.9) / a'-''^ ^ f ^*{a'-''^) ^ {a'-''\^{M)) ^ {<^>*a'^''\M) . 

J<i'(M) JM 

So the pullback is the formal adjoint of the map in this duality pairing. 

Definition 15 (Inclusion map). 5*/ Let A be a subset of A4. The function l : A ^ M defined 
by l{x) = X for every x £ A, is the the inclusion map (or the embedding, or the injection) of A 
into M.. In other words, the inclusion map of a subset of M is the restriction to that subset of 
the identity map on M . 

Definition 16 (Trace operator). /5/ Given two manifolds A4 and M' such that M' C A4, the 
trace operator, tr; 

is the pullback of the inclusion M! ^ M.. If the manifold M is clear from the context, one may 
write tiM' instead of trM,M' ■ If -M' is the boundary of M, dA4, one just writes tr. 

Example 8. Consider an inclusion map l and its associated pullback i* , such as the one depicted 
in Figure^ that generates the inclusion of the manifold A4' , a circle, on the manifold Ai, a disk. 
In local polar coordinates 0' and [9, r) the inclusion map takes the form: 

9 - Lg{9') = 9' 
r = Lr{9') = 1 ■ 
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Figure 9. Pictorial view of the inclusion map t and its associated puUback t*. 



The trace of the 1-form a^^-' — agdO + a^dr is given by: 



tr a' 



(1) dcf^[T6]^^^(i) ^ 



{ae o L){e')d9' 



2.4. Exterior derivative. The exterior derivative, d, plays an important role in differential ge- 
ometry and is defined in the following way: 

Definition 17 (Exterior derivative), \2J^ The exterior derivative on a n-dimensional man- 
ifold Ai is a mapping d : A'^(A^) A'^+^(A^), < k < n — 1, which satisfies: 

(2.10) d(a('=)A6(") =da('='A6(') + (-l)'=a('=) Ad6('), k + l<n, 
and is nilpotent, 

(2.11) dda^*^) := 0, VaC^' e K^{M) . 

For 0-forms, a*-"-* € A'^(A^), the exterior derivative in a local coordinate system {x^, 



jCc"), is 



given by: 
(2.12) 



da(") :^y^dx^ 



For a local coordinate system one has: 



Example 9. In a 3-dimensional Euclidean space and in a local coordinate system (x^,x'^,x^) the 
exterior derivative of a 1-form, a^^"^ — Qida;*, is given by: 

da-" . f - If) d.' A d.= + f - d.» A dx' ■ ^^«' ■ 



dx^ dx^ 



dx^ A da;^ 



\dx^ dx^ J 

And the exterior derivative of a 2-form, fe*^^-* — bidx^ A da;'' + b2dx^ A dx^ + b^dx^ A da;^, is given 
by: 



(2.14) 



d6(2) 



dbi db2 dbs > , , , , ^ 
' da;^ A da:^ A da;^ 



dx^ dx^ dx^ 

Recalling the vector calculus operators, gradient, curl, and divergence, one can see the similarity 



between these and the expressions (2.12), {2. IS) and {2.14), respectively. Indeed there exists a clear 



relation between the two sets of equations, for more details see JIJ \24^ . 

It is possible to show, [1] [76] , that the exterior derivative and the integration of a differential 
form are related by the following result: 

Theorem 1 (Generalized Stokes' Theorem). [76^ Given a k-form a^*^) on a (sub) -manifold 
M. of dimension fc + 1 that is paracompact and has a boundary, then 



(2.15) 



da(^) = / a^'^l 

M JdM 



Two important aspects of this theorem can be stated. One is the fact that it condenses and 
generalizes three key theorems of vector calculus: the fundamental theorem of calculus, Stokes' 
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theorem and Gauss' theorem, for 3-dimensional space, when apphed to 0-forms, 1-forms and 2- 
forms, respectively. The second is the fact that the exterior derivative is the formal adjoint of the 
boundary operator: 

(2.16) d = d^ . 



This property can be obtained by (2.81 and Stokes's Theorem: 

(2.17) (daW,M)<HJ / daW^H^ / <H> (a^, ^A^) . 

Jm JdM 

This is one of the basic identities of the numerical framework to be presented. 

Remark 1. Another important property of the generalized Stokes' Theorem is that if Ai and A4 
are two manifolds of dimension k + 1 with the same boundary, i.e. dAi = dM, we have 

(2.18) / daC^) = /" a^*') = /" a^''^ = [ da^^^ , Va^ e (A'''(7V/() U A'^(7W)) . 
Jm JdM JdM Jm 

Another way of expressing this particular independence of the manifold is that we can always add 

to the manilfold Ai a manifold At with dAA = without effecting any change in the generalized 

Stokes' Theorem. For k — 0, this corresponds to the gradient theorem which states that if two 

points A and B are connected by a curve C, then 

(2.19) j \/(f,-ds = (f){B) - (t){A) . 

This theorem holds for any curve connecting the points A and B , provided that the domain is 
contractible. There is no preferred curve and therefore is seems reasonable to identify all curves C 
which satisfy (2.19). Such an identification will be formalized in the next section. 

The exterior derivative satisfies the following commuting property: 

Proposition 7. Given a mapping ^ : Ai ^ M and the associated pullback, $*, the exterior 
derivative commutes with the pullback: 

(2.20) $*(da('^)) = d$'^(a('^)), Va'*'' e K^{AA). 
and is illustrated as 



1 



A'=(AA) K^{U). 
Proof. For all sub-manifolds ^ in A^, we have 

^ (da«,$(yl)) ^ (a«,9$(yi)) ^ {a^^\^{dA)) = (d<&*aW,y^) . 

Since A was completely arbitrary and it needs to hold for all a'^*'^ € A*'(7W), we have <E>*d = d<i?*. 
See [24l [81] for alternative proofs. □ 

Definition 18. A differential form a^''^ is called exact if there exists a differential form feC^^^) 
such that a^'') — db'^'^^^K The space of exact k-forms is the range of the exterior derivative, i.e. 

B(d; A'-^-i) := dA'"^{M) C A'^(X). 

A differential form a^'^-' is called closed if da^''^ = 0. The space of closed k-forms is the nullspace 
or kernel of the exterior derivative, i.e. 

Z{d:K^) := {Va(*') e h^{M) \ da^''^ = 0} C K''{M). 

It follows from ( 2.11[ ) that all exact differential forms are closed. The reverse is only true on 
contractible manifolds, and is known as Poincare lemma. 

Lemma 1 (Poincare Lemma) . Ifj On a contractible manifold all closed differential forms are 
exact. That is: 

For all a^^) G Z(d; A*'), there exists fo^^'^^) G K^-^{M) such that a^^) = db^''~^^ . 
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Remark 2 (Potentials). The relevance of Poincare's lemma has far-reaching implications. In 
electromagnetics it guarantees the existence of a solution to curlA = B. where B is a given 
magnetic flux density and A is the magnetic potential. It also guarantees the existence of functions 
which give rise to the equations of motion of a mechanical system in Hamiltonian form. Although 
the Poincare lemma is stated in terms of differential forms, the natural isomorphism between 
differential forms and vector fields allows one to state that the Poincare lemma carries over to all 
the existence theorems of potential fields of vector calculus. That is, the existence of solutions for 
all the following equations, for V, A and C : gradV^ — f, curl A = B and divC = /. 



From Definition 17 tlie {n + l)-spaces of differential forms in an n-dimensional manifold A4 



satisfy the following sequence, called the de Rham complex, denoted by (A,d)): 

(2.21) M A"(M) A A^M) A ••• A A"(X) A 0. 

This sequence is exact on contractible domains, in which case S(d; A'""^) = Z(d; A'^). In general 



we have that B{d;A ) C Z(d;A ), as depicted in Figure 10 Hence, one can decompose the 




Figure 10. Pictorial view of de Rham complex on differential forms. The exte- 
rior derivative d maps the elements of A'^ into the kernel of d of A'^+^j B{d; A'^) C 
Z(d;A'=+i). 

space of A;-forms, A''{A4), as: 

(2.22) A'^(X) = B{d; A''-^) ® 6"(d; A'^^^), 

where S^(d; A*^^-'^) is the algebraic complement of B{d; A''^^) in A'^{Ai). One can write any /c-form 
a^''^ as: 

a(fe) ^ dfeC^-i) + c^''\ a^'^) e A''{M), b^'"'^^ G A''-'^{M) and c^'^^ G S"(d; A'=). 

2.5. Hodge-* operator. An important concept for the definition of the Hodge-* operator is the 
the volume form. 

Definition 19 (Standard volume form). fT/ In an n-dimensional manifold A4 with metric gtj, 
a volume form w^"' is defined in a local coordinate system x^ as: 

:= ^\dei{g,j)\dx^ A • • • A dec" , 

where w^") is the standard volume form. 

The Hodge-* operator is then defined in the following way. 

Definition 20 (Hodge-* operator). J71 \24^ The Hodge-k operator in an n-dimensional manifold 
M is an operator, * : A'^(A^) — > A"~'^(7\/J), defined by 

a«A*6W := (a^, , 

where the inner product on the right is defined in Definition \lS\ Application of the Hodge-k to the 
unit 0-form yields *1 := iw^"'. 
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For f,g e C°°{M) and aC'^bC'') e A''{M) the Hodge-* operator satisfies, [SU] 
(2.23a) ★ (/flC') + gfeC^)) = / ★ flC-^) + g * feC^' , 

(2.23b) ★*a('=) = (_i)M«-fe)aW , 

(2.23c) o^^) A = 6^=) A = (flC^), , 

(2.23d) * (a('=) A *5('')) = *(6('=) A ★a^^)) = (a^^), fcC^^) , 

(2.23e) (*aW,*6('^)) = . 

On we have the following relations for the basis forms 

★1 ~ d.TAdy, *dx — dy, -kdy = — dx, -^cdxAdy = 1 . 

It is possible to define an integral inner product of fc-fornis in the following way: 

Definition 21. The space of k-forms, A'^(A^) can be equipped with an inner product, 
(•, •)l2A'=(a^) : ^HM) X A'=(7W) ^ R, given by: 

(2.24) faW,6W) / A . 

Suppose we have a map $ : — > A/". How does the Hodge-T^ transform under this mapping? 
Let us call the transformed Hodge-*, 

Proposition 8 (Transformation of the Hodge-* operator). // a^*^' e K''{M), $ 
A''(A^) and ^ : M —?' Af, then the Hodge operator in M, denoted * : A''{M) -> A"^''(7W), is 
given by: 

★ = $*★($*)-! . 

Proof. The objective is to find a * such that the following diagram commutes, [lO] : 

A'=(A/') — ^ A"-'=(AA) 

1** 1** 
A''{M) — ^ A"-'=(A^). 
Hence * should be such that = Therefore ★ = $* * (<i)*)~^ □ 



The Hodge-* operator enables one to extend the single de Rham complex (2.21) to a double de 
Rham complex connected by the Hodge-* operator 



A"{M) A^{M) ... A A"(X) A 



(2.25) *^ *^ *^ 

0^ A"(X) ^ A"-i(-^) ^ ••• ^ A"(X) M. 
Note the similarity between this diagram and Figure |8] 

Definition 22 (The codifferential operator). 1251 The codifferential operator, d* : A^{jv{) 
A^~^ [M) , is an operator that is formal Hilbert adjoint of the exterior derivative, with respect to 
the inner product {2.24): 

(2.26) (da('-^\b^''A = fa(^-i),d*6W) 

Proposition 9. If a^^-'^^ e A^-'^iM), foC^') e A^{M) and 

(2.27) / d(a('=-i) A *6('^)) = / aC^-^) A *6('^) = , 

JM JdM 

then 

(2.28) d*^^^) := (-l)"('=+i)+i*d*6('^\ V6('^) e A'=(M) . 
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Proof. The proof follows directly from the integrals: 



JM Jm JM 



The first term on the right is zero, (2.27), and the second one is simply: 

'm ^ ' 



M 



Proposition 10 (Nilpotency of codifferential). The codifferential satisfies 



□ 



□ 



Proof. It follows directly from (|2^, ( |2.23bp and ( [2lT |. 

Due to the Proposition [lOj the (n + l)-spaces of differential forms satisfy the following sequence: 

(2.29) i— AO(X) ^ A\M) ^ ■■■ ^ A"(X) ^ R. 

Definition 23. A differential form a^*^^ is called co-exact if there exists a differential form feC^+i) 
such that a^*^' = d*^'''"^"'^'' . The space of co-exact k-forms is given by 

B{d*;A''+^) ■.^d*A''+\M) cA''{M). 

A differential form a*^*"'' is called coclosed if d*a'^^^ = The space of coclosed k-forms is given by 

Z{d*;A'') := {VaC') e A'^(X) | d*a^''^ = 0} C A\M) . 

Proposition [l0[miplies that S(d*; A*-') C Z(d*; A^^-i). Only when 6(d*; A'=) = Z(d*; A'^^^) is 
the sequence in ( 2.29[ ) exact. One can decompose the space of fc-forms, A'^(A^), as 




jB(d*:A'--') 
1 2'(d*;A*) 



2:(d*;A'=) 



Figure 1 1 . Pictorial view of the complex of differential forms acted upon by the 
codifferential operator. The codifferential d* maps the elements of A*^ into the 
kernel of d* of A'^-^, B{d*-A'') C Z(d*; A'^'-i). 

A'^iM) = Bid*;A''+^)(BB''{d*;A''+^), 

where B''{d*; A'^+^) is the algebraic complement of ;B(d*; A'^^^) in A''{A4). And one can write any 
fc-form a*^''^ as: 

^{k) ^d*6('^+i) +0^=), flW e A'^(X), feC^'+i) e A''+\M) and c^'^) G B^id* ; A''+^). 

Definition 24 (Laplace-DeRham operator). fJl \TE^ The Laplace-de Rham operator, A : 
A'^(A^) — > A'^(A^), is a map given by 

(2.30) Aa^^') := (d*d + dd*)a('=), Va^^' G A'=(M) . 

Proposition 11 (Self-adjointness of Laplace-de Rham operator). Under the conditions of 



(2.21) the Laplace-de Rham operator satisfies 

(Aa«,6('=) 



= faW,A6«" 



i.e. A* = A. The Laplace-de Rham operator is self adjoint. 
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Proof. The proof follows directly frora the definition of the Laplace-de Rham operator (2.30) and 
( p^ . □ 

2.6. Hodge decomposition. A fc-form is called closed if da^' '^ = and exact if there exists a 
{k — l)-forni foC^'^i) such that a^'^-' = d6''''~^-'. Since do d = 0, ( |2.11[ ), every exact form is closed. 
For manifolds which are non-contractible the Poincare Lemma [1] states that not all closed forms 
are exact. 

One can define the space of harmonic forms as those differential forms which are closed, but 
not exact, 

(2.31) •H''' :=Z(d;A'^)nS"(d;A'=-i) =^ Z{d; A'') ^ n''{M) U B{d; A''-^) . 

For contractible manifolds H'^ — 0. Similar to (2.22) we can decompose A''{Ai) into its nullspace 
and its complement, 

A''{M) = Z(d; A'=) ® Z"(d; A*^). 
If we combine this decomposition with (2.31), we obtain the following Hodge decomposition, 

(2.32) A''{M) = B{d; A''-^) ®n'' ® Z"(d; A*^). 
Proposition 12. If ty a^^'^ = or dM = 0, then 

(2.33) Z{d*;A'') = B%d;A''-^) Z=(d; A'=) = S(d*; A'=+i) . 
Proof For all a^^^ € Z(d*;A'=) and all b'^''-^^ e A''{M) we have 



(aW,d6('=-i) 



L2A'=-i(X) 

Therefore, ± S(d;A'=-i), i.e. a^^) S S^(d;A'''-i) = 6'=(d; A'^-^). Therefore, Z(d*;A'=) C 
^c(j.^fc-i)^ Conversely, if a^^') G ^-^(d; A'^'-^), then (d*a('=\ fe(''"^))^2A.(^) = for all b'-''-^^ 
which implies that d*a('=) = 0, therefore S=(d;A'=-i) C Z(d*;A'=). So we have B^id; A''-^) = 
Z(d*;A'=). 

For aU a^^) e Z(d;A'=) and ah b^''^ e B{d* ; A''+^), i.e. there exists a c^^+i) such that fot'^) = 
d*c('=+i). Then 



L2A_fe+i(A^) 

This implies that Z'=(d; A'=) = ^(d*; A'^+i) 



(a(^),d*c('=+i)) 



□ 



Using Proposition 12 in (2.31) yields 

H'^' = Z(d;A'=)nZ(d*;A'=), 

or 

(2.34) •H'= = {VaGA'=(X) |da = 0, d*a = 0}. 

Harmonic forms are therefore both closed and coclosed. Using Proposition 12 and (2.32) allows 
us to write the Hodge decomposition as 

(2.35) A'=(X) =S(d;A'=-i)®-H'=®S(d*;A'^'+i). 

Remark 3. Note that the Hodge decomposition \2. 32 ) is true, whether Ai has a boundary or not, 
whereas {2.35) is only true if dA4 = 0. The Hodge decomposition in the form (2.32) will play an 
important role in the remainder of this paper. 

Corollary 2 (Hodge decomposition). P]/ Let M. be a compact boundaryless oriented Rie- 
mannian manifold. Every e^'^-' G A'^{M) can be written in terms of a^^^^^ G A'^^^{M), fof'^+i) g 
A^+'^{M) and c^''^ G H'' such that 

(2.36) e^^) =da('=-i)+d*6('=+i)+cW . 

Remark 4. Let M be a compact boundaryless oriented Riemannian manifold and a^^^ G H^ . 
Then it follows that c''^-' is the harmonic solution of the Laplace-De Rham operator, 

(2.37) ( d*c^''^ = and dc'^^^^ = ) ^ Ac^^'> = 0. 
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Remark 5. Although da'^''-^\ d*6('=+i) and e'^^^ are unique, the differential forms a^'^ and 
^(fe+i) ^j,g generally not unique, because we can replace a'^'^"^-' by a'^^~^^ + dp^'^"^-* for any p^*-'"^) 
and b'^'^^^^ by feC^+i) + f\* q{^+'^) Qi^d this will also satisfy the Hodge decomposition. 

Remark 6. The k-th de Rham cohomology group of M, , is defined as: 

, ^ Z(d;A^) 
■ 6(d;A'=-i)' 

It is possible to prove that the space of harmonic k-forms, Ti,^ , is isomorphic to the cohomology 
group H'' , see Moreover, the de Rham theorem, see J75| /, states that for a finite dimensional 
compact manifold this group is isomorphic to the homology group defined in algebraic topology, 
the isomorphism being given by integration. This connection between differential geometry and 
algebraic topology plays an essential role in the development of the numerical scheme presented 
here since it enables the representation of differential geometric structures as finite dimensional 
algebraic topological structures suitable for using in the discretization process, as it will be seen in 
the following sections. 

For a manifold with boundary, the Hodge decomposition theory takes a somewhat different form. 
A fc-form a'*^^ G K''{M.) is called parallel or tangent to a manifold Q C if tr^_g(*a'^'^)) = 
and is called perpendicular or normal to a manifold Q C if tr_A4,g(a'^''''') — 0. In this way the 
following spaces can be introduced 

(2.38) ^t{M) = {a^^) e K^{M) \ a^^^ is tangent to dM} , 

(2.39) A^(7W) = e A'=(7W) I a^^) is perpendicular to dM} . 

Note that for the space of harmonic forms da'^'"'' = and d*a'^'^^ = is a stronger condition than 
Aa^'^'' = when M has a boundary. 

Proposition 13 (Hodge decomposition theorem for manifolds with boundary). 'J / Let 

M be a compact oriented manifold with boundary. The following decomposition holds 

(2.40) A''{M)^dA';;-\M)®n''{M)®d*A^+\M). 

2.7. Hilbert spaces. On an oriented Riemannian manifold, we can define Hilbert spaces for 
differential forms. Let all fi{xi^, . . . ,Xi^) be functions in L'^{Ai), then a^*^-* be a fc-form in the 
Hilbert space L^A'^(r2) is given by 

a^'"^ = ^ fi{xi^,. . . ,Xi,J dxi^ A dxj2 A ... A dx^^. 

i 

The norm corresponding to the space L'^k^{M) is ||a''°^ ||^2A^fc = {^'^^\ '^^^'') j^k (^j^y Although 
extension to higher Sobolev spaces are possible, we focus here to the Hilbert space corresponding 
to the exterior derivative. The Hilbert space HK^{M.) is defined by 

and the norm corresponding to HK^{M) is defined by 

Ik^'-'llk'^ :=||aW||i.^. + ||da«|li.^.+,. 

The de Rham complex, or the Hilbert version of the de Rham complex, is the sequence of maps 
and spaces given by 

E HA°{M) A HA\M) ^ i?A"(X) 0. 

An important inequality in stability analysis, relating both norms, is Poincare inequality. 

Lemma 2 (Poincare inequality). Consider the de Rham complex (A, d), then the exte- 
rior derivative is a bounded bijection from Z'^(d, L^A*^) to S(d, L^A'^'^^), and hence, by Banach's 
bounded inverse theorem, there exists a constant cp such that 

(2.41) ||a||^A. <cp||da|L2A., Va G Z^(d, L^A'^), 
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Which we refer to as the Poincare inequality. We remark that the condition B{d, i^A*^ ^) is closed 
is not only sufficient, but also necessary to obtain this result. 
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3. An introduction to Algebraic Topology 

An important step in the our numerical framework is the discretization of a manifold (physical 
space) in which the physical laws are embedded. By this we mean the partitioning of a manifold 
in a collection of non-overlapping subspaces (cells) such that their union is the whole manifold 
under study. This partitioning into a set of distinct subspaces yields a representation of a manifold 
in terms of a finite number of points (0-cells), lines (1-cells), surfaces (2-cells), volumes (3-cells) 
and their analogues of dimension k (/c-cells) . A collection of fc-cells of such a partition is called a 
A;-chain. Several types of geometrical objects (fc-cells) can be used for the subdivision. In this work 
we will focus on quadrilaterals and their generalizations to higher dimensions, (singular fc-cubes). 

In this section we intend to formally introduce the representation of these collections of k- 
cells and to associate them with a suitable algebraic structure that enables a correct discrete 
representation of the original manifold. The branch of mathematics that provides such a formal 
discrete representation of a compact manifold is algebraic topology. 

Having defined a discrete representation of the manifold in terms of fc-chains, we can assign 
values to the various space elements by defining the dual space of the chains space; the so-called 
k-cochains. Duality pairing between chains and cochains allows us to introduce operations on 
cochains as formal adjoint of operations on chains. In this way, the introduction of operations on 
fc-cochains mimics the operations on /c-forms discussed in the previous section. 

The main reason this introduction on algebraic topology is given here, is that we will explicitly 
employ the close relationship between differential geometry in a continuous setting and algebraic 
topology in the discrete setting in the following sections. This relationship is also employed by 
many others, for instance, [SI [501 HZl [ZH]- For a more thorough treatment of algebraic topology 
the reader is referred to [ ^ Hi l H5 1 [5^ [75] . 

In this paper we assume that the dimension of all manifolds is finite and dim(A^) = n. 

3.1. Cell complexes and the boundary operator. 

Definition 25 (fc-cell). \28\. l5l| / A k-cell, ^k): of a manifold M of dimension n > k is a set of 
points of A4 that is homeomorphic to a closed k-ball Bk = {x g M'' : ||x|| < 1}. 

The boundary of a k-cell T(^k): d'''(k)} is the subset of Ai associated by the above mentioned 
homeomorphism to the boundary dBk — {x E M.'' : ||x|| = 1} of Bk- 

The topological description of a manifold can be done in terms of simplices, see for instance 
[52l [75l l82] or in terms of singular k-cubes, see [44l [45l |80j. From a topological point of view 
both descriptions are equivalent, see |19j . Despite this equivalence of simplicial complexes and 
cubical complexes, the reconstruction maps to be discussed in Section [4] differ significantly. For 
mimetic methods based on simplices see [H [U [T51 [3T1 [S51 [SS] , whereas for mimetic methods based 
on singular cubes see [?! 1^ IT71 [M[ [B0[ [71] . 

Here we list the terminology to set a homology theory in terms of /c-cubes as given by [IS] : 

M = real line. 

/ = closed interval [—1, 1]. 

M''' = M X M X • • • X M (fc factors, k > 0) Euclidean /c-space. 
l''^=lxlx---xl{k factors, fc > 0) unit fc-cube. 

By definition l'^ is a space consisting of a single point. 

Definition 26 (Singular /c-cube). A singular fc-cube in a n-dimensional manifold Ai is a 
continuous map T(^k) ■ — ^ Ai, < k < n. 

Definition 27 (Degenerate singular fc-cube). '4-5'J If for the map Tf^f.^ : /'^ — > Ai, < k < n 

there exists ani, 1 <i < k such that t^j.) (xi, X2, . . . ,Xk) does not depend on Xi, then the singular 
k-cube is degenerate. 

Corollary 3. Any non-degenerate k-cube is a k-cell as defined by Definition \25\ 
Remark 7. A non-singular k-cube is a submanifold of Ai according to Definition^ 
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Remark 8. Select an orientation in M'^' , then the determinant of the Jacobian of the maps 
determines the (inner) orientation of the k-cell according to Definition^ 



- cell 




Figure 12. Example of a 0-cell, a 1-cell, a 2-cell and a 3-cell. 



Figure 12 depicts some examples of /c-cells in a manifold M — R . The fc-cells are geometric 
objects which represent the geometric objects shown in Figure |8] Before we formally define the 
boundary of /c-cubes, we first define faces of a fc-cube. 

Definition 28 (The faces of a singular fc-cube). f^5[ [761 / For < k < nO let T(j,) be a singular 
k-cube in Ai. For i = 1,2, . . . ,k, we define the singular [k— I) -cubes ^iT(fe_i) 
M , by the formulae (face maps ) 

AiT(^k~i){xi,X2, ■ ■ ■ ,Xk-i) = T(fe)(xi, . . .,Xi^i, -l,Xi, . . . ,a;fc_i) 

BiT(k^l){xi,X2, . . . , Xk-l) = T(k){xi, . . . ,Xi^l,+l,Xi, . . . ,Xfe„l) 

AiT(/j_i) is called the front i-face and i?iT(j,_i) is called the back i-face of T(^i^y 



Figure 13 depicts some examples of faces of /c-cells in a manifold Ai — . 






Figure 13. Examples of the faces (in dark) of a 1-cell, a 2-cell and a 3-cell in . 



Remark 9. Note that the face maps, A^rj^.i) and -BiT(fc_i) defined in Definition 28 are inclusion 
maps as defined in Definitional^ 

Definition 29 (The boundary of a singular fc-cube). 145^ The boundary d of a singular 
k-cube T(^k)y k > is given by 

k 

(3.1) 



This definition describes the boundary of the submanifold r(fe), see Corollary [TJ 

Definition 30 (Cell complex). f2R A cell complex, D, in a compact manifold M is a finite 
collection of cells such that: 

(1) (D) is a covering of M.. 

(2) Every face of a cell of D is in D. 

(3) The intersection of any two k-cells, t^^i^'^ and ai^i^-^ in D is either 

• T(i;) and (T(fe) share a common face; 



T(k) 



n cr 

n a( 



(fe) 



or: 



T(k) I I O (fc) 
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In this framework we will focus on: points (0-cells), curves (1-cells), surfaces (2-cells) and k- 
dimensional generalizations (fc-cells), considered as non-degenerate fc-cubes. Figure [m] depicts 
an example of a cell complex in a compact manifold C M^. The above presented definitions 
constitute a formalization of the concept of discretization of space. 



Manifold 



Cell Complex 




Figure 14. Example of a cell complex. Left: a three dimensional compact man- 
ifold. Right: the k-ce\\s that constitute the cell complex. 



Definition 31 (fc-ciiain). !28^ Given a cell complex D, the space of k- chains of D, Ck{D), is 
the free Abelian group written additivively, (see 145}/ ), generated by a basis consisting of all the 
oriented, non- degenerate k-cells of D. A k-chain C(fe) in D is an element of Ck{D). 

A fc-chain, C(j,) G Ck{D), is a formal sum of /c-cells, T(fc).i G D: 

Vc(i,) e Ck{D) C(fe) = ^ cV(fc),i, r(fc),j e D . 

i 

Formally, given a set of fc-cells, it is possible to generate any fc-chain by specifying the coefficients 
in the chain. Although this is possible for arbitrary fields, we will, in the description of geometry, 
restrict ourselves mainly to chains with coefficients in Z/3 = {—1,0,1}. The meaning of these 
coefficients is : 1 if the cell is in the chain with the same orientation as its orientation in the cell 
complex, -1 if the cell is in the chain with the opposite orientation to the orientation in the cell 
complex and if the cell is not part of the chain. The orientation of the cell is implied from the 
orientation of M*^ and the map as pointed out in Remark [s] For an example of fc-chains, see 
Figure [15] 



Cl = + ''"(1),4 + TCl),!! + T(1),12- 

-T"(l),6 - T"(l),3 - T"(l),8 - T"(l),7 



f-2 = ^(2).} + T(2).2 + •'"(2),3 +'^l2).4 



T-(()),i T(o,_,i 




Figure 15. Example of a cell complex, a 1-chain and a 2-chain. 

Definition 32 (CoefHcient vector for chains). The space of k-chains, Ck(D), can be rep- 
resented by a column vector containing only the coefficients of the chain. That is, there is an 
isomorphism tp: 

(3.2) ib:CkiD)^RP, p = rank(Cfc(i^)) , 

defined by 



(3.3) 



V'(c(fc)) = V'(2jcV(fc)^,) = [c^ • --cPf, p = rank{Ck{D)) , 



MIMETIC FRAMEWORK 



27 



where the rank of Ck{D) is the number ofk-cells in the cell complex D and the c' are the coefficients 
of the k-chain C(fc-). The k-chain, C(-j,-), is printed in boldface, whereas the vector C(^i;^ of coefficients 
is printed in regular face. 

We can now extend the boundary operator applied to a fc-cell, Definition |29[ to the boundary 
of a fc-chain. 

Definition 33 (Boundary of a fc-chain). I2}A The boundary operator, d : Ck{D) — > 
Ck-i{D), is an homomorphism, defined by 

(3.4) 9c(fc) = 9^ cV(fc)^i := ^ {T(k),i) , 



where the action of the boundary operator on the k-cubes is given in Definition 29 



The boundary of a fc-cell T(fe) will then be a (fc — l)-chain formed by the faces of r^^). The 
coefficients of this (fc — l)-chain associated to each of the faces is given by the orientations. 

i 

with 

= 1, if T(j,_x),i has the same orientation as the face of T(j.) ^ ; 
e'j — —1, if T(fe_i) i has the opposite orientation as the face of t^j.) j ; 
— 0, if T(fc_i) i is not a face of T(fc) ^ . 
And the boundary of a 0-cell is 0. 

(3.5) 9c(fc) = a(^c'T(fc)j) = ^c?dT(k-),j = ^c''e}r(fc_i),, . 



Example 10. The boundary of the 2-cell T(^2).i i"^ Figure 15 has the following boundary: 

One clearly sees that, in this case, we have: e\ — 1, e\ = 1, e\ = —\, el = —1 and e\ = for 
i^ {1,2,7,9}. 

Definition 34 (Incidence matrix for ciiains). The coefficients e* constitute an rank(Cfe_i) xrank(Cfc) 
incidence matrix E(k~i,k) with {E(k-i,k))ij = ej. 

Proposition 14 (Boundary operator and incidence matrix). Let ?/;(c(j,-)) = C(fc), then 

V'(5C(fc)) = E(^k-l,k)C(k)- 

Proof. 

'/'(oc(fc)) W ■,/; 2_^cJejT(fc_i)_, = 2^(^ej = E(^k-i.k)C(k) ■ 

\ ij J ij 

□ 

An essential result that follows from the definition of the boundary operator is that applying 
the boundary operator twice results in a zero chain. 

Theorem 2 (The boundary of the boundary is empty). The boundary operator satisfies: 
(3.6) ddc^k) = 0(fc-2), Vc(fc) e CkiD). 

Proof. The proof follows from repeated application of Definition [29] to fc-cells which extends by 



Definition 33 to fc-chains. See also [151 [75] for this derivation. □ 



This simple result will have profound consequences for the discrete operators to be defined. 



9 o 9 = is illustrated in Figure 16 
Proposition 15. 



E(fc-i,fe)E(fc_fc_i_i) 
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A 



da 



dda 



Figure 16. The boundary of the boundary of 2-cell is zero because all edges 
have opposite orientations. 

Proof. For aU C(fc_|_i) e Ck+i{D) we have 

„Th. [2],.„„ , Prop. [14] , .Prop. [14] , , 



□ 



15 



-(2) 



'^(2).! 



r{2),2 



Example 11. Apply the boundary operator to the 2-chain in Figure 

T(2),3 + '''(2),4- 

9C(2) = 9(r(2)_i + T(2)_2 + T(2),3 + 'r(2),4) 

= + '^(1),4 + '^(1),11 + T(l),12 - '^(1)^6 - Til),3 - Til).,8 - Til).,7 ■ 

This is the 1-chain depicted in the middle in Figure [75[ Applying the boundary operator again we 
get: 

ddC(2) = d{T(i)^i + T(i),4 + T(i)^ii + T(i)42 " '''(I), 6 ^ '''(1),3 " ■''(1),8 " = • 

This result could also be obtained by the use of incidence matrices, since they are a matrix repre- 
sentation of the topological boundary operator: ^{k-2,k-i)^(k-i.k) — 0- Consider the cell complex 
in Figure 15 The incidence matrices ^(o,i) and ^(1^2) are 



-(0,1) 



E(l,2) - 



1 

-1 





-1 



1 









1 

-1 






-1 



1 






which produces E(-o4)E(x.2) — 0, as expected. 

Definition 35 (Cycles and boundaries). fSUf A k-chain C(j.) for which dc(^k) = 0(*:- 
a cycle. A k-chain Cf^^) which is the boundary of a (fc + l)-chain b(j,_)_2), i.e. C(-j,-) = 
called a boundary. 



is called 



dh 



(fc+i) 



IS 



The space of fc-cycles in a cell complex D is denoted by Zk{D) and the space of ^-boundaries 
in D is denoted by Bk{D). Theorem [2] implies that every boundary is a cycle, Bk{D) C Zk{D), 
but the converse is generally not true. Therefore, consider the factor space Hk{D) consisting of 
those cylces which are not boundaries 



(3.7) 



Hk{D) 



Bk{D) 



Hk{D) is an equivalence class and any two cycles C(j,) and d^f,) are associated with the same 
element whenever the difference C(j,) — d(j.) is a boundary. H^i^D) is called the homology group 
and two elements which differ by a boundary are called homologous, [75]. In many applications. 
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dim(iJfe(Z?)) is finite and this dimension is called the k-th Betti number. Two spaces can only be 
topologically the same, if the their respective Betti numbers are the same. 

Example 12. Consider the cell-complex depicted in Figure [T7| which contains a 'hole' in the 
middle. In FigureU^all 0-cells have positive orientation by default. The orientation of the 1- and 



\ 





/ 




7 8 
3 4 




/ 





\ 



Figure 17. A non contractible cell complex. 



2-cells are indicated in the figure. The incidence matrices which relates the 0-cells to the 1-cells 
and the 1-cells to the 2-cells, are given by 



-(0,1) 








-1 
1 











-1 



1 





E(l,2) = 



The matrix ^(q^i) describes the connectivity between the 8 nodes and the 12 line segments. The 
range of the matrix is spanned by its 12 column vectors. Not all these column vectors are linearly 
independent. The rank of this matrix is 7, hence the range of this matrix has dimension 7. Since 
the dimension of the null space of matrix ^(1.2) is 8, there is one element in the null space o/E(i 2) 
that is not in the range o/E(o.i). This confirms that on a non- contractible domain, not every cycle 
is a boundary. The dimension of the homology group is in this example equal to 8-7=1. This 
corresponds to the number of 'holes' in the domain (formally known as the Betti number). The 



-1) 



harmonic chain corresponding to Figure 17 is given by 

h(i) = (1,-1,0,0,1,1,-1,1,-1,0,0, , ^ 

The topological Hodge decomposition of the space of k-chains is given by 

Ck = Bk®Hk® Zl Zl = Ck\Zk . 

The harmonic k-chain that belongs to the space Hk is defined as 

h(fe) = {9h(fc) = I Ja(fc+i) e Cfc+i, such that h(fc) = da(^k+i)}- 

Remark 10. Note that can only be obtained from global considerations. Therefore, the de- 
composition Ck = Bk ® Hk © Z^ can only be obtained globally. 

The boundary operator defines a differential graded algebra of degree -1, {Ck{D),d), from all 
the /c-chain spaces of the cell complex, giving rise to a chain complex. 



Definition 36 (Graded algebra of degree -1, {Ck{D),d) ). I28f A chain complex is a family 
{Cfc(D),9} of k-chains and boundary operators, that constitutes a sequence 



(3.8) 
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Only when all H''{D) = is this sequence an exact sequence. 

In the previous section, we discussed how differential forms and operators acting on these 
differential forms behave under a mapping Similar constructions are also feasible in algebraic 
topology. These maps are called chain maps and cochain maps. 

Definition 37 (Cliain maps). Let $ : A4 ^ JV be a homeomorphism between manifolds A4 and 
M , then this map induces a homomorphism between k-cubes in Ai and k-cubes in M given by 




where <I>jj (T(fc)) = <i> o r(s.) for any singular k-cube T(^k) '■ l'^ ^ -M, for k — 0,1,2, ... ,n. $jj 
is extended to chains to be the homomorphism $(j : Ck{D^) — > Ck{Dj\f), the induced chain 
homomorphism, by putting 



{k),i) ■ 

It can be shown, [U], that $jj maps non-degenerate fc-cubes in M into non-degenerate fc-cubes 
in A/" and that the diagram 



I *4 



Ck-i{DM) < Ck{DM)- 

commutes. This follows directly from the fact that the face maps. Definition [28j commute with 
the chain map: $tt [At*") = Ai ($jtr'=) and $(j (BiT'=) = Bi ($ttr'=). Therefore, we have that 

(3.9) ao$j=$joa, 

the chain map $jj commutes with the boundary operator. The image of the boundary is the 
boundary of the image. This is the discrete equivalent of Proposition [2j 

Since the boundary operator, d, commutes with the chain map, $[t maps the cycles Zk{DM) 
into Zk{Dj\/) and the boundaries Bk{D_\4) into Bk{Dj^), for all fc = 0, . . . , n. Therefore, a chain 
map does not alter the topology. In particular, the incidence matrices for the cell complex Dm 
and Dj^ are identical. 

3.2. Cochains. 

Definition 38 (Cochains). J^] The space of k-cochains, C^{D), is the space dual to the space of 
k-chains, Ck{D), defined as the set of all the linear maps (linear functionals), c'^^^ : Ck{D) — >■ M, 
and we write 

(3.10) (cW,C(,)) :=cW(c(,)), 
to represent the duality pairing. 



Remark 11. Note the resemblance between [3.10) and (2.8). This similarity forms the basis of 
the mimetic framework discussed in this paper. 

Remark 12. Strictly speaking k-cochains are homomorphisms from Ck{D) into M.: C^{D) := 
Hom{Ck{D);m.). 

Remark 13. Since the cochain spaces are the dual spaces of the chain spaces, we might as well 
write -.^Cl. 
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Proposition 16 (Canonical basis fc-cochains). Given a basis of Ck{D), {r^^^yj} with i — 
1, • • • ,p with p = rank[Ck{D)), there is a dual basis of C^{D), {r^'"''*}, with i = I,-- - ,p with 
p — rank{Ck{D)), such that: 

(3.11) rW.Xr(fe),,)=<5i-. 

All linear functionals can be represented as a linear combination of the basis elements: 

i 

The proof of this proposition can be found in any book on algebraic topology, for instance, |52j . 
With the duality relation between chains and cochains, we can define the formal adjoint of the 
boundary operator which constitutes a sequence on the spaces of fc-cochains in the cell complex. 

Definition 39 (Coboundary operator). '28'J The coboundary operator, 6 : C*'{D) — > (7*^+^(1?), 
is defined as the formal adjoint of the boundary operator: 

(3.12) (<5c('=),C(,+i)) := (c('^),ac(,.+i)), VcC') e C\D) and Vc^k+i) e Ck+i{D) . 
Proposition 17. The coboundary operator is nilpotent, 

(3.13) Mc^^) ^ 0, VcC^) e C^{D) . 

Proof. The proof follows directly from Theorem [2j Vc^*") e C^{D) and Vc(fe) e Ck{D) 

□ 

Definition 40 (Cocliain complex). JW^ A cochain complex is a family {C^{D),5) of k-cochains 
and coboundary operators, that constitutes a sequence 

(3.14) ^ C''-\D) A C''{D) A C''+\D) A • • • . 

Analogous to Definition |32] we have 

Definition 41 (Coefficient vector for cochains). The space of k-cochains, C^{D), can be 
represented by a column vector containing only the coefficients of the chain. That is, there is an 
isomorphism tp: 

(3.15) V' : C'iD) MP, p = rank(Cfc(L')) , 
defined by 

(3.16) Vi(c('^) = V^(Ec»^^'^'') = [ci---Cp]^, p^rankiC'^iD)) , 

i 

where the rank of C^{D) is the number of basis k-c ochains in the cell complex D and the Ci are 
the coefficients of the vector c^^^ in Proposition 16 The k-cochain, c'^^\ is printed in boldface, 
whereas the vector c*-*^-* of coefficients is printed in regular face. 

Proposition 18 (D uality pairing in terms of coefficients). Duality pairing of k-cochains 
with k-chains, (3.10), in terms of the coeff dents, 'ij}[c'^'^^) = c^'^-* and ?/'(c(j,')) = C(fc), is given by 

(3.17) (cW,C(,.)) = (cW)'^C(,). 
Proof. 

□ 
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Proposition 19 (Incidence matrix for coboundary operator). Let ^(c^*^^) — then 

(3.18) = (E(,,fe+i))^ . 

Proof. For all c^*^) G C^(D) and for all C(^k+i) e Ck+i(D) we have 

m 



(E(fc,fc+i)C(fe+i)) 



(E(fe,fe+i)) 



C(k+1) 



k) (k)\ 



□ 



Proposition 20. Let he the incidence matrices for the coboundary operator, then 

(3.19) ^(k+i,k)^(k,k-i) _ 

Proof. 



£(fe+l,fc)^(fc,fc-l) Prop- [191 



-{k,k+l 



(fc-l,fe) j 



-{k-l.k)^{k,k+l) 



T Prop. M5^ 



0. 



□ 



Recalling the duality pairing of differential forms and manifolds (2.8) and the duality pairing 



of chains and cochains (3.101 one clearly notices the algebraic similarity between the continuous 



and the discrete worlds. Expression (3.121 is nothing but a discrete Stokes' theorem. 



Definition 42 (Cocycles, fc-coboundaries and cohomologous cochains). f7g| / The cochains 
cC^) for which Sc^'^'^ — 0*^*^+^^ are called cocycles. The set of all k-cocycles is denoted by Z^{D). 
A k-cochain that can be written as the coboundary of a {k — l)-cochain, c^'^-' = 5d^''~^\ is called 
a /c-coboundary. The space of all k-coboundaries is denoted by B^{D). 



Corollary 4 (Cochain space decomposition). From Proposition 11 is follows that B^{D) C 
Z'^{D). We therefore consider the cohomology group 

Z^{D) 



(3.20) 



of all k-cocycles which are not k-coboundaries. The spac e H^ {D) is the topological equivalent 
of the space of harmonic forms T-0'{M.) defined in Section 2.6 The space of k-cochains can he 
decomposed as 



Compare this decomposition with the Hodge decomposition given in (2.32) 



Proposition 21 (Cochain well-posedness). Consider the following cochain relation, da.'^^^ — 
f(/c+i)^ yjllli g^{k) g {^Z^Y 0.''^^ e B^^^ . Then there exists a solution aS^^ and the solution is 

determined up to a cochain in Z^ . 

Proof. There exists a solution a'^'"'\ because fC'+i) g 5^^+^. Now, suppose there exists two solu- 



tions, a^''^ and a^^\ such that 5a^^^ = f^^'+i) and (Ja^"^ = fC^+i). Then 5{a^^^ - a^"') = oC^+i) 



(k) 



Jk)-, 



Since a 



(fc) 



a^''^ e Z'', the solution for feC^) = fC^+i' is unique for a^'^) S {Z''^. 



□ 



Example 13 (Example 12 continued). The dimension of the cohomology is equal to the dimen- 
sion of the homology due to the duality pairing. The harmonic cochain corresponding to Figure 1 7 
is given by 

h^i) = a(l, -1,0, 0,1,1,-1,1,-1,0,0,-1)^, a e M. 
This harmonic cochain models a circulation around the hole in Figure with strength a. 

Remark 14. Note that the harmonic cochain is global. Only global considerations lead to the 
determination of the harmonic cochains, which is the main reason these solutions cannot be rep- 
resented by local methods, such as finite difference, finite volume or finite element methods. 
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Definition 43 (Induced cochain map). Let $jj : Ck{D_\4) Ck{Dj^), k = 0,1,..., n, then 
for every k-cochain. ci'^'' € [D j^f) , there exists a k-cochain Cm' G C''{Dm) such that 

(c«,$„C(fe)) = (cW,C(fe)) := ($»cW ,C(fe)) , Vc(fe) e Cfc(i?M) , 
tiiii/i f/ie cochain map defined by Cm"* = $'cl''''. 

Remark 15. The cochain map satisfies (cm\^'uC(fe)) = ($'*cl'^\ C(j.)) . Compare this relation 
with the duality between a space map and its pullback as discussed in Proposition^ 

Corollary 5. The cochain map commutes with the coboundary operator S, 

(3.21) ^o$» = $»o(5. 



<J)» 



Proof. 

(5$«cW,c(,)) °'=5.^($»cW,ac(,)) °'=^'3^(cW,$,ac(,)) <P (cw,a$,c(,)) 



□ 



Remark 16. The commutative property So^^ = ^^oS is the discrete analogue ofdo^* = $*od, 
see Proposition^ 

Since the cochain map commutes with the coboundary operator, maps cocycles onto 
cocycles and fc-coboundaries onto fc-coboundaris, for all fc = 0, . . . , n. This property ensures that 
the discrete Hodge decomposition retains its structure when we deform the cell complex with a 
chain map. 

3.3. Dual complex. With every cell complex D one can associate a compatible dual cell-complex, 
D, [79]. Before we can define a dual cell-complex, the adjoint of the faces should be defined. 

Definition 44 (Faces and cofaces). J7P| / The faces of a k-cell are those [k — l)-cells that form 
the boundary of the k-cell. The cofaces of a k-cell are those (k -f \)-cells which have the k-cell as 
common face. 

Example 14. Consider two neighboring rooms in a hotel. The two neighboring rooms have the 
wall in between them as common face. Therefore, the cofaces of the wall between the rooms are 
the two rooms. 

Definition 45 (Dual cell). With every k-cell, ^k); *^ ^ there corresponds a [n — k)-cell, f(^n-k)7 
in the dual complex D, such that the dual cells of the boundary of a k-cell are the cofaces of the 
dual cell T(„-/c) • 

Definition 46 (Dual cell-complex). A dual cell-complex, D, is the smallest cell complex ac- 
cording to Definition 30. that contains all dual cells f^n-k) o.s defined in Definition 45 



If we denote the map, which associates a cell r^j.) with its dual T(„_fe), by *, i.e. *T(^k) = '^(n-fe) 



and if we denote the coface operator by d* then Definition 45 states that the following diagram 
commutes. 



Therefore we have *d — d**. An example of these relations is shown in Figure [18) 
We make a distinction in cell-complexes with and without boundary. 
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Figure 18. Pictorial view of Definition 45 The dual of a line segment in D C 
is gray shaded surface in D. The cofaces of this gray surface are the two adjoining 
volumes in D. The boundary of the line segment are the two endpoints in D. The 
dual of the two endpoints are again the volumes in D surrounding these points. 



Definition 47 (Dual cell complex without boundary). // the collection of all dual cells of 
D, denoted by *D , constitutes a cell complex according to Definition \3C\ then the dual cell complex 
D = *D and both the cell complex D and its dual D are cell complexes without boundary. 

Corollary 6. For a cell complex without boundary the cell-complex D is dual to the cell complex 
D and the cell complex D is dual to the cell complex D modulo orientation. As a consequence 
exists and is = ±*. This allows us to express the coface operator in terms of the dual operator 
* and the boundary operator d by 



(3.22) 



* d = d** 



d* — *9* = ± * (9 
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Figure 19. Example of a primal grid and dual grid is shown (not necessarily a cell 
complex), with for every cell Tj-^) the associated dual cell f(„_j,) and corresponding 
orientation. On the right the staggered grid, an overlap of both grids, is shown. 



Remark 17. The dual grid in Figure [7ff| is only a cell complex in case the top side is connected 
to the bottom side and the left side to the right side. This makes this cell complex on the surface 
of a torus. 

Remark 18. Note that the primal cell complex, D, is chosen outer oriented. Then by duality, 
the dual cell complex, D, is inner oriented. In fact, the orientation itself does not change, only 
the corresponding cell changes. This was shown before in Figure^ 

Remark 19. // we equip the cell complex, D, with an outer orientation, the the dual complex, 
D, models geometric objects with inner orientation. Alternatively, inner orientation could be 
represented on the cell complex D, in which case the outer representation is modeled on the dual cell 
complex D . In this respect, dual cell complexes are able to model the inner- and outer- orientation 
as discussed in Section W^ 
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3.4. The boundary of cell complexes. A collection of fc-cells forms a cell complex if the 
boundary of these fc-cells are also in the cell complex. Therefore, cell complexes consist of an 
interior Di and a boundary part Df,, where D — Di U Di,. The boundary part contains cells up 
to degree (n — 1). Since D is dual D, we have that Di U Dh is dual to Di U D/,. More precise, 
Di = *Di and Di, — *Di,. The individual parts do not need to be cell complexes themselves. In 
case Di is a cell complex on a manifold with boundary, then Di is not a cell complex, because 
not all faces of the fc-cells in Di, k = 1, ... ,n, are also in Di. Examples can be found in [7S] and 



Examples 15 and [T6| below. 



Definition 48. (The boundary of cell complexes) Let D be a cell complex. If *D is not a 
cell complex, let D be the smallest cell complex which contains all the k-cells of *D. All the k-cells 
in D \ *D form a {n — 1) -dimensional cell complex, Df, := dD called the boundary of D. The 
dual cells of dD with respect to the (n— 1)- dimensional embedding space, Dh := dD = *{dD) form 
a {n — \)- dimensional subcomplex of D, called the boundary of D. A dual cell in Df, is given by 
ff^n-i-k) = The boundary and its dual are a cell complexes, because the boundary of the 

boundary is empty and Definition \4'l\ 

Example 15 (ID primal- and dual- cell complex). Let the interior part of the primal cell 



complex, Di, be given by three 0-cells, connected by two 1-cells, see Figure 20 Then its dual, 
Di, consists of two 0-cells and three 1-cells. The outer two 1-cells are open ended, and therefore 
Di is not a cell complex. To make it a cell complex, Di must be closed by adding the boundary 
cells, Db. The boundary part Di, consists of two 0-cells. The dual of Db are the boundary cells of 
the primal cell complex. In this one- dimensional example, Db also consists of two 0-cells. Note 
that for the primal cell complex, the boundary cells in Db coincide with the cells in Di. Although 
the boundary points in Db were already contained in Di , the addition of orientation requires us to 
treat the boundary points as distinct points from Di . This is because outer orientation of a point 
in a ID embedding space differs from outer orientation of a point in a OD embedding space, see 
Example^ So Di and Db together form the primal cell complex D, and Di and Db together form 
the dual cell complex D. 



D =^ =9!^ =^ ^ >_ ' — ■ > ' D 

Figure 20. Primal and dual cell complexes, split into their interior and boundary 
part. The corresponding orientation to all fc-cells are indicated. Note that the 
end-points in D have two types of orientation. 



Example 16 (2D primal- and dual cell complex). Now consider in two dimensions the 
interior part of the primal cell complex, Di, as given in Figure \2l\ The corresponding dual is 
the interior part Di, of the dual cell complex D. Then D becomes a cell complex by adding the 
boundary cells Db. Note that there do not exist 0-cells at the corners. The dual of Db are the 
boundary cells Db. They again coincide with the cells in Di. Although Db is already contained 
in Di as points and lines, the role of the points and line segments is completely different. This 
is indicated by the orientations in Figure \21\ The outer orientation along the boundary - a ID 
cell complex - differs from the outer orientation of these same geometric objects considered as 
elements from Di embedded in 2D. 

Remark 20. The boundaries Db and Db are boundaryless cell complexes. 

The construction of the boundary of a cell complex and its dual given above, indicates that 
the boundary of a cell complex is disjoint from the cell complex itself. In order to formally 'glue' 
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Figure 21. Interior and boundary parts of the primal and dual cell complexes, 
D = DiU Db and D = DiU Db, respectively. Orientation of all cells are included. 



the boundary to its cell complex, we introduce the chain inclusions ij : {dD) — {D) and 
ijj : dCk{D) Ck{D), such that, for = 0, . . . , n — 1, 

(3.23) (r(fc) e Ck [dD)) = t^^) e Ck [D] , and ^ {ff^k) e Cu{dD)) = G Cfc(^) . 

Chain inclusions «j and are the discrete analogues of the continuous inclusion map defined in 
Definition |15[ Using duality pairing between chains and cochains, we can define the associated 
cochain maps : (D) — >• {dD) and : C''{D) C^{db) of the inclusion maps. 

Definition 49 (Discrete trace operator). For all c*^*^-' e C^{D), there must exist € 
C^{dD), such that 



(3.24) 

The map is then defined by 
(3.25) 



c''"':«tj(c(fc))) , Vc(fc) e Cfe (9D) 



The cochain map is defined similarly. The cochain maps ^^ a nd are the discrete analogues of 
the of the pullback of the inclusion map given in Definition 16 We therefore will write tr and tr 
instead of andiK 

Remark 21. Note that the inclusion chain maps can be defined more generally: Let K be a 
sub-cell complex in D, then we can define the map ij : C]~{K) — > Ck{D) and its associated dual 
operation on cochains. Although useful in some applications, in this paper we restrict ourselves to 
K = dD. 

If we associate the primal complex D with a cell complex endowed with outer orientation and 
the dual complex D with a cell complex endowed with inner orientation, we have 

Definition 50. (Tangent k-cochains) A k-cochain c^'^^ e C^{D) is caHerf parallel or tangent 
to the cell complex D, if ti (c^'^-') = 0^'^'-'. We will denote the set of all tangent k-cochains on a 
given cell complex D by C^{D). 

Definition 51. (Normal k-cochains) A k-cochain c^'^^ € C^{D) is called perpendicular or 

normal to the cell complex D, if ii (c^*^-*) = 0^'^-'. We will denote the set of all normal k-cochains 
on the cell complex D by C,^ [D) . 
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If the dual cells are labeled with the same number as the associated primal cells, and if the 
orientations of primal and dual cells are in agreement with the orientation of the embedding space, 
then the incidence matrix E(p_]^ p) on the dual complex is related to the incidence matrices on the 
primal complex by 

(3-26) ^{p-i,p) = E.f^_p^n_p+i), p=l,...,n. 

Remark 22. This seemingly uninteresting relation has quite some consequences for numerical 
operators to be developed. The incidence matrices encode the boundary operator, d, but by Stokes 
equation, also the coboundary operator, 5. The coboundary operator is the discrete analogue oj the 
exterior derivative. Take for instance p — \ in {3.26), then E^^''^^ = ^-^ is the discrete analogue 
of the exterior derivative acting on 0-forms. We identified this operator with the gradient operator. 

(3.27) E("--i) = (E(„_i,„))^ ^ E(o,r) = (e(1'0) ' 

The coboundary operator acting on {n — l)-cochains discretely represents the exterior derivative 



acting on (n — l)-forms, i.e. the divergence operator. Therefore, (3.21), states that the discrete 



divergence operator (on the primal complex) is the transpose of the discrete gradient (on the dual 
complex). In vector calculus it says div = — grad . The minus sign is a consequence of the fact 
that orientation is not taken consistently into account in vector calculus. Furthermore, vector 
calculus does not reveal that these two operators act on spaces with a different type of orientation 
(inner and outer). The notion of inner- and outer- oriented complexes will lead to staggered grids. 

Remark 23. Consider in 2D the discrete gradient operator, applied to t he d ual cell-complex. Then 
the incidence matrix E.^^'*^^ relates the 0- and 1-cells as shown in Figure 22 Note that the discrete 



gradient does not give rise to 1-cells between the 0-cells on the boundary. 



Figure 22. A sub-cell-complex showing all 0- and 1-cells that are involved in 
the discrete gradient on the dual cell-complex in 2D. It consists of four internal 
0-cells (•), eight boundary 0-cells (■) and twelve line segments. 
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4. Mimetic operators 

The general philosophy in this section is to develop projections and coprojections, denoted by 
TT and TT*, respectively, onto a finite dimensional space of differential forms, A^(0;Cfc), such that 
an operation T at the finite dimensional space Af^{ft;Ck), behaves in the same way as in the 
continuous level, A'^(Sl). By 'behaves in the same way' we mean that either tt o T — T o tt or 
t:* o T — T o TT* ^ i.e. when one of the following diagrams commute 



Af^in) 



A'(n) 

-1 



In the previous sections, we denoted a manifold with either Ai or M. In computational engineering 
it is customary to denote the manifold - or computational domain - by fJ. Therefore, from now 
on we will refer to the manifold as fl. 



4.1. Reduction, reconstruction and projection operator. A discretization involves a pro- 
jection operator, tt, from the complete space A''{il) to a subspace A^(ri;Cfc) C A'^(i7). In this 
subspace we are able to express differential forms in terms of fc-cochains defined on fc-chains, and 
corresponding fc-cochain interpolation forms (often called basis functions or basis forms). There 
is a clear relation between the continuous representation using differential geometry and the dis- 
crete expressions in algebraic topology. There exists a reduction operator, TZ, that integrates the 
fc-forms on fc-chains to get fc-cochains. On the other hand, there is a reconstruction operator, X, 
to reconstruct fc-forms from fc-cochains using appropriate basis forms. These mimetic operators 
were already introduced before in the work of Hyman and Scovel |33j and in Bochev and Hyman 
[S]. A composition of the two operators gives the projection operator tt = X o TZ as illustrated 
below. 

A'=(r!) A^(r!;Cfe) 



n 



These three operators together set up the mimetic framework. The coprojections cannot be 
written directly in terms of reduction and reconstruction. However coprojections are expressed 
in terms of Hodge-* operators and projections. The latter are written in terms of reduction and 
reconstruction. In this section we focus on the discretization of the differential forms, involving 
the reduction, reconstruction and projection operators. 

Definition 52. The reduction operator TZ : A'^(f2) C^{D) is a homomorphism that maps 
differential forms to cochains. This linear map is also called the De Rham map and is defined by 
integration as 




(4.1) 



{TZd 



(k) 



Then for all C(j,) G Ck{D), the reduction of the k-form, ' G A (17), to the k-cochain, ' € 
C^{D), is given by 

(4.2) aW(c(,)) {TZa^'^\c,,,) ^^^EU ^ e'(7^aW, r(,),) ^ V, 



,(fe) 



,(fc) 



It is the integration of a fc-form over all fc-cells in a fc-chain that results in a fc-cochain. Note 
that since T(^k),i is compact and C^{D) is a finite dimensional free Abelian group, the reduction 
operator TZ is well-defined for a'^^^ = a{xi, . . . , a;„)dx*i A ... A da;*'' with a(xi, . . . , Xn) G L\oc{^)- 
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Definition 53 (Integration over the domain ft). Let a'") be any n-form, then the integral 

f a(") := (7^a("),a;(„)), 
where the chain <jJ{n) = '^iT{n).i (so all c* = +1) covers the entire computational domain Q. 



We can list the following properties of the reduction map. First of all, the reduction operator 
is a non-injective and surjective map. 

Definition 54. Consider two k-forms a^^\a^2^ e A'^'(57) and a^'"' ^ a2^\ they are in the same 
equivalence class if 



0. 



(4.3) 7^al"^=7^4'^ ^ 7^(4'^-af) 
Secondly, the reduction map commutes with respect to differentiation. 

Lemma 3. The reduction map has a commuting property with respect to continuous and discrete 
differentiation, 

(4.4) 7^d = (57^ onA''{n). 
This commutation can be illustrated as 



Proof. This property can be proven using Stokes' Theorem (2.15) and the duality property (3.12), 

(4.5) (7^da(^),C(,)) e f daC^) ^ f a^') O {na^-\dc,,,) ^ (SUa^'Kc,,,) . 

□ 

Thirdly, the reduction map commutes with the pullback, $*, and the cochain map, 

Lemma 4. Let $ : flj^ flj^ be a continuous map between two manifolds, let $j : Ck{Dj^) — > 
Ck{Dj\f) be the associated chain map of k-chains between two cell complexes and let Dm be o, 
covering of the manifold fl^ o,nd Dj\f the covering ofQj^/, according to Definition \3(^ Then the 
reduction map commutes with the continuous and discrete pullback, 

(4.6) 7^<I>* = $«7^ on h^{VLj^). 
This commutation is illustrated by 

K'^iSiM) ^^{^M) 

where <I>* is defined by Definition 14 and is defined in Definition 43 
Proof For aU C(fc) e CkiDj^) and for all a^^' e A''{flj^) we have 

□ 

The operator acting in opposite direction to the reduction operator is the reconstruction oper- 
ator, I, which maps fc-cochains onto finite dimensional fc-forms, and is defined as follows. 
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Definition 55. The reconstruction operator I : C^{D) K^{Q.;Ck), also called the Whitney 
map, is an isomorphism that maps cochains back to differential forms. The reconstructed differ- 
ential forms belong to the space Af^{^l^,Ck), which is a proper subset of the complete k-form space 



K^{Vl). While the reduction step is clearly defined in Definition 52, in the choice of interpolation 
forms there exists some freedom. Although the choice of a reconstruction method allows for some 
freedom, I must satisfy the following properties: 

• Reconstruction X must be the right inverse of TZ, so it returns identity (consistency prop- 
erty), 

(4.7) ni = Id onC^{D). 

• Like TZ, also the reconstruction operator T has to possess a commuting property with respect 
to differentiation. A properly chosen reconstruction operator X must satisfy a commuting 
property with respect to the exterior derivative and coboundary operator, 

(4.8) dX = Xd onC'^iD). 
This commutation can be illustrated as 



• The reduction should commute with continuous and discrete pullback, i.e. 
(4.9) 1$' = on C''{Dj^). 

This commutation relation can be represented by 



• Let Di = *Di be as defined in Definition 48 For all c'" '^^ € C" ^{Di) there should exist 
a c'^^^ e C^{Di), such that 

(4.10) ic("-'=) = Zc('=\ 

where X reconstructs cochains on Di and X reconstructs cochains on the cell complex Di . 
The same must hold in the opposite direction. 

Both should also hold for the boundary part of the pr imal and dual cell complexes, i.e. 



forDb andDb- LetDb = be as defined Definition 48 For all c^"^^^''^ e C"~^ '"(-Db), 
there should exist a c^*^) £ C^{Db), such that 

where X reconstructs cochains on D^ andX reconstructs cochains on D^. Again the same 
must hold in opposite direction. 

Remark 24. Moreover, we want X to be an approximate left inverse of TZ, so the result is close 
to identity (approximation property), 

(4.11) XTZ = Id + 0{hP) inA'=(17), 

where 0{hP) indicates a truncation error in terms of a measure of the grid size, h, and a polynomial 
order p. 

The composition in the last remark gives rise to the definition of a new operator, tt/j, which is 
a projection operator. 
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Definition 56. Define the operator tt^ ■ A'^(il) Af^{n;Ck) as the composition I o TZ. It allows 
for an approximate continuous representation of a k-form a^''\ 

(4.12) 4^'' = ^^aC-^) = X7^a('=^ nha^^^ e Al{n) C A'^ift). 

where ITZa^^^ is expressed as a combination of k-cochains and interpolating k-forms. 

Proposition 22. The operator iTh : A'^{il) — ^ Af^{n; Ck) is a projection operator. 

Proof. For tt/j to be a projection, it must be a homomorphism, which is true since both the 
reduction and the reconstruction operators are homomorphisms, therefore for all a^''\ 6^'^-' £ A'^(J7) 
it holds 

(4.13) 7r,(aW+6W)=7r,aW+^,6W, 

and the projection operator must be idempotcnt. Let aj^'"'' e A^(il;C'^), then i^hOh'^ — o.h'\ 
and so tth = Id on A'^{n;Ck). For aU a|f^ G A{;(fJ;Cfc), there exists aC") e A'=(f7) such that 
a\^'> = TT/iflC^) ^ina^''\ so 

(4.14) I7^4'=) = X7^(I7^a('=)) <^ X7^a('=) = a|f\ 

□ 

Definition 57 (Bounded projection). For a projection operator tt^ to he a useful operator, we 
require it to be a hounded projection operator, i.e. for Ci < oo, 

fA TK\ II II hha'^^'^Wf," , ^ 
(4-15) \\T^h\\c(KKAr)-= sup II (fc)ii <Ci, 

where || • Ha*" ^•s some norm defined on A'^'(57). 

Corollary 7. Alternatively, one can write for the bounded projection operator 

(4.16) IIWIIa- <Ci||aW||A.. 

h 

From the triangle inequality applied to ||a(''') || a.^ it follows that there exists C2 = Ci + 1 < 00, such 
that 

(4.17) ||(/rf-7r,)a(^)|U. <C2|laW|U.. 

h 

Proposition 23. Using the projection operator, the differential form space A''{n) can be decom- 
posed into a projected space and its complement, so A''" — Ajj ® ^h '^ , o.i^d so any k-form can be 
uniquely decomposed as 

(4.18) a(^) -4''+4''^'" = ^^a('^) + (W-^^)aW, ^a^'^^ e A\n). 

Note that (Id — Tih) is also a projection, hut now onto the unresolved part. As a consequence of 
the direct sum decomposition and Proposition 22. iTh o (Id — tth) = [Id — tt/j) o tt^ = 0. 

The projection is not an orthogonal projection and therefore the complement space i s n ot 
orthogonal to the projected space, so AJ^^'*^ ^ A^'"^. This will be demonstrated in Example 
the next section. As a consequence the projection, tt^, is not self-adjoint, 

(4.19) (^/.a('=\M'=)) ^ {a^'\n>,b^'^^) , ya^'^\b^''^ e A'=(17). 

Both a^'^) e A'^(J7) and its projected part iTha'^^^ G A^Jil; Ck) are in the same equivalence class, 
as defined in Definition [54] so 

(4.20) TlTTha^'"'^ ^Tla^'"'^ and 7^(/d - 7r,,)a('=) 0, 

and for the special case of integration of a volume form over the whole domain, this integral 
preserving property of the projection gives 
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As for the reduction operator, also the projectfon operator is non-injective and surjective. Then 
because of Definitions [54] and [56] there also exists an equivalence class for the projection operator. 
The projection of a'*^' is not unique, but depends among others on the underlying cell complex, 
as is illustrated in the following example: 



Example 17. Consider a uniform and a Gauss-Lobatto grid, then 7r""'a^ 7r| ' e A^([— 1, 1]) , 
but 7r|J'"a('^^ 7^ n^a^^\ see Figure 23 for an example of a 0-form and a 1-form. The difference 
between the two projections reduces with grid refinement, since 



(fc)i 



<l(/- 



]4.11| l 



0(hP). 



Moreover note that, according to Proposition 

(fe) ^ ' 



TT^'aC^) and so 7r™'7r,f 
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it can be observed that Ti'^J'^TT^'a'-'^'' = 




(a) zero-form (b) one-form 



Figure 23. These figures illustrate the non- uniqueness of the projection operator 
for zero-forms (a) and one- forms (b). In (a) the dots indicate the 0-cochains, in 
(b) the line segments are the 1-cochains. The resolutions are intentionally kept 
low to show the differences. 



In Section]3]staggered cell complexes were introduced. The two complexes are dual with respect 
to each other as defined in Definition ]46] and illustrated in Figures [20] and [21] This allows us to 
define two projection operators, one projecting onto the subspace A^ft; Ck) and one projecting 
onto the subspace A'^(f7;C'fe), where Ck = Ck{D) and Ck — Ck{D). 

Definition 58 (Canonical projection operators). Given the outer- oriented cell complex D 
and its dual inner- oriented cell complex D, we define two projections: 

• TT/i = XTZ, where the reduction, TZ, and the reconstruction, I, are performed on the interior 
part of the primal cell complex, i.e. Di. The corresponding subspace A^(r2; Ck) is the space 
of finite dimensional k-forms, with k-cochains associated with the outer -oriented cells in 
Di, such that 

A^(r!;Cfe) := 7r,,A^-(^^) - I7^A'=(^!). 

• TT/i = TTZ, where the reduction, TZ, and the reconstruction, X, are performed on the interior 
part of the dual cell complex, i.e. Di. The space A^j(r2; Ck) is the space of discrete k-forms, 
with k-cochains associated with the iimei- oriented cells in Di, such that 

Atin;Ck) 7^ftA^-(f^) -i7tA'=(r!). 

To every projection we can define a corresponding coprojection. 
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Definition 59 (Canonical coprojection operators). Given projection tt^ using the outer- 
oriented cell complex Di and projection tt/j using the inner- oriented cells in Di, we define two 
coprojections: 

• TT^ is defined as = (— l)'^'^"~*''^*7r/i*. This coprojection is defined in terms of a projection 
of {n — k) -forms, ★a*^'^\ on the inner- oriented (n — k) -chains in Di and then represented 
in terms of the k-form basis functions on the outer- oriented cell complex Di. The corre- 
sponding subspace is TT*^A''{il). 

• TT^ is defined as tt^ = (— l)'^("~*'')*7r;i*. This coprojection is defined in terms of a projection 
of {n — k)-forms, -ka''^\ on the outer- oriented [n — k)-chains in Di and then represented 
in terms of the k-form basis functions on the interior part of the inner- oriented dual cell 
complex, i.e. Di. The corresponding subspace is 7r^A'^(f2). 

Remark 25. A coprojection is also a projection according to Proposition \2S\ 

Remark 26. Although projections and coprojections in Definitions \58\ and \5S\ were defined with 
the interior parts of cell complexes D and D. The same projections and coprojections can be 
defined for the boundary parts of the cell complexes D and D, i.e. Db and D}y. 

The projection and coprojection operators possess commutation relations with the operators to 
be defined in the next subsection. The projected differential forms differ, because the forms are 
reduced on different chains as illustrated in Example 23 See also Figure 31 in the next section. 

Proposition 24 (Equivalence spaces of differentials forms tt/iA*^ and tt^A''). Define the 
spaces 

nnA' {ai'^ \3a^''^ G A\n) s.t. ai^'^ = Tr.at^') } , 

and 

^*A'= := {4''^ |3a('=) G h^{n) s.t. a[^^ = nl^aC'^ } , 
then K\ = iihA^ = ^^h^^' ''^ general iiha^^^ ^ T^h^^'^^- 

Proof . The reducti on o f t^to*-'^-* on (ti— fc)-chains in Di yields a (rt— fc)-cochain in Di and according to 



(4.10) in Definition 55 there exists a fc-cochain on the cell complex D such that Tc^" ^'^ = Xc^^\ 



Ifaf eA^ Pro^ 



therefore 7r,*A C n^A . By a similar argument one can show that n^A C tt^A , so we have 
TT/iA'' = Tr^A*". That TTha^'''' ^h^-^''^ will be shown in Figure 31 in Section [Hj □ 

Proposition 25 (Equivalence projections on subspaces). If a'^'^'^ G Ajj(i7;Cfc) then Tr^a'^'^-' — 
Tiha^^'^ and if a^^^ G Al{n;Ck) then n^a'^^'^ = n^a^''^ 

Proof. 

a\^^ = T^hO-^h ^ because tt/j = Id on Aj^ = tt/jA'^ 

= <a'-J^^ because tt* ^ Id on A^ = tt^A'^ 

Therefore ajj'^'' ~ T^hO-'ii^ — '^h^'ii^ ■ The proof on the dual complex is the same. □ 

Remark 27. Although tt^o''^-' = Tr/jc''^-' holds in the finite dimensional subspace A'^(fl;Ck), it 
does not holds on the entire space A'^(ri). The difference between these two projections could give 
rise to natural and derived operators, f^jj. However, from our point of view projections tt^ and 
coprojections tt,* are essentially different operators. 

Definition 60. Whenever we refer to the projection n it can be any of the projections from 
Definition\58[ 



Although TT/i = Id on A^, (Proposition 22 1, it will be shown that it can be a very useful 
operation when changing the expression of a finite dimensional fc-form. A finite dimensional k- 
form a"j^^ is usually expressed in terms of a A:-cochain a^'^^ G C^{D) and interpolating A:-forms. If 
this is not the case, but the A;-form a\ ' is expressed in terms of an Z-cochain corresponding to a 
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chain in Q :— Ci{D), where D is either the primal ceU complex D or the dual cell complex D, 
then a projection is used to express the finite dimensional fc-form in terms of its corresponding 
fc-cochains. For this special case we introduce a separate projection operator, ttm- 

Definition 61. Define a special, bijective projection ttm : A^(i7; Ci) — )■ A^j(f2; Ck), such that 

(4.21) TTM ^ Id onKl{n;Ci). 

So TiMO'h = cLh, but its expression changes in terms of cochains, from an l-cochain in C\D) to a 
k-cochain in C''{D), and its expression changes with respect to the basis functions. In terms of 
reduction and reconstruction this is 

(4.22) Oh = ina = ITZiilZa) = iruah- 

The seemingly redundant projection 7rj\/ will reappear almost everywhere, but we will not ex- 
plicitly mention this in this section. It will be more instructive to show its action with concrete 
reconstruction operators in Section [Sj 

4.2. Discrete operators. In this section we discuss operations of the operators discussed in 
Section[2j restricted to the set of finite dimensional subspaces A^, i.e. d, d*, Ah and (•, ■)h- 



4.2.1. The exterior derivative. First consider the exterior derivative. With (4.4) and (4.81 a com- 
muting property of the projection with respect to the exterior derivative can be shown. 

Lemma 5. The projections tt/j = TTZ and tt^ — TTZ commute with the exterior derivative, 

(4.23) dTT/i — TT/jd and dnh — ifhd on 

This can be illustrated as 

1' 



Proof. Express the projection in terms of the reduction and reconstruction operator, then for all 



dTT/ifl = dTTZa ISTZa ITZda = Tr^da. 
The proof for tt/i = ITZ is the same. □ 

4.2.2. The Hodge-k operator. The projection and coprojection possess a commuting diagram prop- 
erty with the Hodge-Tk-. 

Proposition 26. We have for all a'-'^'' G A'^'(57) the following commutation relations, 

ir'^-k = -kiTh and ★ tt^ = -k^-k . 

Proof. From Definition [58] we have 

TT^ — (-1)'^'("~''') ji-f^i, <^=> TT** — -kn^ and * tt^ = tt/^ ★ . 

So the Hodge-Tk- operator at the continuous level commutes with the Hodge-* in the finite dimen- 
sional setting with respect to the projections tt/j and tt^. 

k * ^ A n~k 



A" — A 



^k 








K - 





□ 
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Remark 28. Similar commutation relations can be set up on the dual cell complex and this gives: 
for alla^''^ G A''{n) 

with the associated commutating diagrams. 



Corollary 8. Proposition 26 combined with Proposition 25 gives that for a^''^ e 



Prop.g6^ ^ Prop. [25| 



'★TT 



h 



In Definition 20 the Hodge-* was defined as a mapping from fc-forms to {n — /c)-forms. As a 



consequence, not only the dimension of the chain over which they are integrated changes, but also 
the orientation changes from inner to outer orientation or vice versa. This corresponds to the 
vertical relations in Figure [8j On finite dimensional subspa ces, the Hodge-* operator is a mapping 



58 



from Af^{Q; Ck) to ''(il; C„„fc), as defined in Definition 

With the Hodge-* and exterior derivative, the finite dimensional double De Rham complex can 
be set up similar to ( 2.25| ) as 



■X 







4.2.3. The codifferential. With Definition |58| a commuting property of the coprojcction with re- 
spect to the codifferential can be shown. 

Proposition 27 (Commutation coprojections with codifferential). For all a'^^^ e A'^(f2) 
we have 



d*7r;* and Tr^d* 



Proof. 



(-1) 



n+k + l 



The proof for vr,* is the same. 



□ 



This means that the codifferential is exact with respect to the coprojection, i.e. the following 
diagram commutes 



Afc+i 



A^^ 



Aj+i 



A'= 



d* 



Corollary 9. Combining Proposition 27 and Proposition 25 shows that for a^^^ £ A^(f2) ive have 



^ Prop. [25] ^ *Prop. |27]^ ^ Prop. [25]^ 

4.2.4. r/ie pullback. The projection commutes with the puUback as follows. 



Lemma 6. Let it be either Wh or ith as in Definition 60 For all a^^^ € A''{n_\f) and for $ 



Hm ^ ^J\f ! there exists a commuting property between the projection operator tt and the pullback 
such that 

(4.24) $*7r = 7r$* on A'^iflj^). 

This commutation can be illustrated as 



A^i^M) 
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Proof. The proof is based on the commuting properties of the discrete puUback, : 



An ahernative proof can be given, based on (2.9) and (4.20). Integrate over a k-ce\l r^j,), then 



$*7ra 



(fe) 



,(fc) 



l |4.20[ 



□ 



This commutation does not hold for the coprojections. Combining the commutation relations 
( |4.23[ ) and ( |4.24[ ) gives 

7r<J>*d = $*7rd = $*d7r = d4>*7r = d7r$* = 7rd$*. 
4.2.5. The wedge product. The next operator to be considered in the subspaces C A*^, < 



k < n, is the wedge product. A product of a k- and ?-form from subspaces Af^ and A^ gives a 
(k + Z)-form that is not in the subspace A^"*"'. It therefore requires again a projection step. 

Definition 62. A discrete wedge product is introduced .such that Ah ■ A^ x A'^ — > A^^', given by 

(4.25) a«A.6i'):^vr(ai'^)A5i')), 

where tt is either tt^ or tt/j , as defined in Definition \6C^ 

As a consequence the discrete wedge product, A^, approximates the wedge product. A, because 



Ahb^ = {Id-7:)iarAb^h') 



(fc) 



0{hP 



Let us verify that the discrete wedge product satisfies the same properties as the original wedge 
product. Let a'"i^\c\^^ € A^, bH"^ € A'^ , with 2k + I < n. Using the linearity of the projection it is 
straightforward to show that 



'^h Of^ - Ah Oft + Ah Oft . 



Also the skew-symmetry follows from the linearity of the projection. 



„(fc) A Ai) _ 



- A b^) = (-1)-. {b^^ A a^) - i-lM' Aft . 



The third property of the wedge product is associativity (2.2b). As stated already in [18], the 
associativity property is in general not satisfied. Now let a e AJJ, 6 G A'^, c G A™, with k + l + m < 
n, ther[^ 

(a Aft 6) Aft c- a Aft (6Aftc) 

= TT (7r(a A 6) A c) - TT (a A Tr{b A c)) 

= TT [{a A b) A c - {I - TT){a A b) A c ~ a A {b A c) + a A {I - -K){b A c)] 
a A {I ~ TT){b A c) - {I - 7T){a A b) Acl = 0{hP). 



0{hP) 



0{hP) 



When considering bilinear products, the normal wedge product and discrete wedge product are in 
the same equivalence class, 



(4.26) 

In case fc + / = n, we get 



7^ 



(4'=)Aft6|^))=7^(4'=)A6«). 



ai^^Ab^\ 



^Both sub- and superscripts are intentionally suppressed for readability. 
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The discrete wedge product in combination with the exterior derivative gives the Leibniz rule 
( [2I0I ), 

d(a« A, = dvr(a('=) A l^^) = .d(aW A 

(4.27) =.(da,f A5« + (-l)Mf'Ad6i'' 

Since the puUback operator $* commutes with the projection, the discrete algebra homomorphism 
is satisfied 

4.2.6. T/ie discrete inner product. Last operator to be defined is the inner product restricted to 
the finite dimensional subspace Ajj. 

Definition 63. Define a discrete inner product (•,-)^ : x Af^ AJJ as 

(4.28) {a['\€\--n{{a^^\b[^^).^-^}, 

with either vr = tt/j or tt = tt/j . T/iis discrete inner product is bilinear, symmetric, positive definite. 
Corollary 10 (Discrete inner product). The corresponding L2 inner product {■r)L'^iih ■ 
A^ X Ajj ^ M is essentially the same as {2.24), because of {4.20) and that Aj; C A*^. Let a^y{\b^y{^ € 
A^, then 



(4.29) 



(o<"i.l.'>).<"> = («<">,(,!.'>)_ 



The Hodge-* was defined in Definition [20] as a combination of the inner product and the wedge 
product. Now when using the discrete wedge product and discrete inner product, the Hodge-* 
remains unchanged, and therefore 

(4.30) {ai'\b^,^\ = a\^^A,.bi^\ 



In weak formulations, integrals over are considered. In that case (4.30) reduces, according to 
(|4.26|) and (|4.29|, to 



(4.31) (ai'=),5i^-))^ 



4^)a*5|^-). 



This result shows that in a weak formulation, the original inner product, wedge product and 
Hodge-* operator can be used. 

4.3. Discrete Hodge decomposition. Let A^ C A'^ be the space of finite dimensional differ- 
erential forms, and let (A/i,d) be a finite-dimensional subcomplex with dA^^ C A{^^^. Since the 
exterior derivative commutes with the projection operator, it follows that ;B(d,A^) C S(d, A'^), 
Z(d,AfJ C Z(d, A*"') and that S(d,AfJ C Z(d, A^). We can make the following decomposition of 
the space of discrete differential forms, 

A^ = S(d,A|;-i)®6^(d,A^i). 

The space of discrete harmonic forms is defined as = A/'(d, A^) n S'^(d, Aj'j^^). This gives 
the following decomposition, S^(d,A^-^) = V.'^ © 2:'=(d,A^). Although 6^(d, A^) C 6=(d, A*''), 
in general, "H^ (;t H'' and Z'=(d,A^) ^ Z=(d,A'=). They both depend on d*, since Z'=(d,A'=) = 
S(d*,A''+^), of which we know that it does not commute with the projection tt/i on A*^, see 
Remark [sj Both harmonic spaces 'H'^ and Tif^ are finite dimensional spaces and their dimension 
depends on the topology of the domain fl. Because the dimension only depends on the topology 
(its Bettie number), we have dim — dim Ti''. The gap between and H'' vanishes as /i — >■ 0. 
See Theorem 3.5 in |4j for the definition and details about the gap between the two harmonic form 
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spaces. Substituting the decomposition into the previous decomposition gives the discrete Hodge 
decompostion, 

= S(d,A|;-i)®Hf,eS(d*,A^+i). 

Example 18. Consider the Hodge decomposition, Corollary^ of a k-form, a^^^ € , in terms 
o/fe(fc-i) e A'^-i, /iC^) G andc^''+^'> £ A^+\ 

aW =d6('=-i)+/i(^)+d*c('=+i). 

Then apply the projection operator, tt^, 

^,a(^) = d^,fe(^-i) + n.ih^'^^ + d*c(^+i)) = d6|f-i) + (4'=) + d*ei'+\ 

ee(d,A^) ee<=(d,A^) 

where /if G Uf^ ^ and e|f+'^ G KI+\ hut e^'+'^ ^ ^^cC^+D. 
//we apply the coprojection operator, tt^, to a*-*^-*, we obtain 



ee<=(d*,A;;) eB(d*,A;;) 

where h)^J G H^. a«rf r)^""'^ G K^-', but r)^''' ^ tt*1 

The continuous Hodge decomposition, discrete Hodge decomposition and the cochain space 
decompostion, are related to each other by means of the reduction and reconstructino operators. 

Proposition 28. From the definition of the reduction operator TZ, the Hodge decomposition for 
k-forms and the cochain space decomposition, it follows that 

(4.32) 6(d; A'=-i) A B^ n'' A i^^ Z'=(d;A''0 A {Z'')". 

From the definition of the reconstruction operator X, the Hodge decomposition for finite dimen- 
sional k-forms and the cochain space decomposition, it follows that 

(4.33) B'' A B{d; Afr^), iJ^' A nt {Z'^f A Z%A- Afj. 

A consequence of this proposition is discrete well-posedness. 

Theorem 3 (Discrete vk^ell-posedness). For g S(d;A'') and a^^^ G Z'=(d;A''), the 

equation da'-'"'-' = /('^+^) is well-posed in the sense that there exists a solution and the solution is 
unique. Then the discrete equation da,f — fl''^'^^ is also well-posed if and only if the projection 
TT is either tt/i or TTh- 

Proof. For /C^+i) g B{d;A'^) and a^''^ G Z'^(d;A'^) there exists a unique solution by the Hodge 
decomposition (2.32 1. The discrete equation follows by projection, 

= 7r(da('=) - /('=+!)) - dTra^'^) - Tr/t^+D - d^") - /f 
Moreover we can write 

^daW - TT/C^'+i) - I(<57^aW - 7^/('=+l)) = ^(JaC^) - f^^+i)). 



By Proposition [28) 7^a('=) = a^^) G {Z'')" and 7^/('^■+l) = ff'^+i) G B''+^. Because I is bijective, 
the discrete problem is also well-posed, since the cochain relation (Ja*^*^' = jg well-posed, see 

Proposition 21 Furthermore, it follows that a^^"^ has a unique solution, because a*-*^^ is a unique 
cochain of the problem fe^'"'^ — fC^+i). □ 



Example 18 showed the decomposition of the finite dimensional fc-form a\^ ^ G A^. The following 
example explains how to extract the harmonic part from this /c-form. 
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Example 19 (Discrete harmonic forms). Following Example\l^ let Tr/ifl^'^^ € and h^^^'' e 
T-L^ its harmonic part by the discrete Hodge decomposition. In order to find the harmonic form 

(k) ' ' 

hf^ , we first determine the harmonic chains, ii(k) the cell complex, as described in Example 



12 



and the harmonic cochains, h^''\ as described in Example IS Both the harmonic chains and 
cochains can be obtained by purely topological considerations and the calculations are performed 
through the incidence matrices. For dim{'H'^) = 1, we set h^^^^ ~ I {ah.^^^) , where a is obtained 
from the definition of reduction, 

(4.34) /7^a«,h(,A.(.h«,h.-A 



(,); = ^.n ,n(,,/ _ - ^^^^^^^^^^^ 
When dim{H^) > 1, we repeat this process for all harmonic chain-cochain pairs and set 

dim(«;:) 

(4.35) /.(f) = Yl 

The harmonic form h^'^} can be found in a similar way as in the case with tt/j. Find the {n ~ k)- 
harmonic chain-cochain pairs, (h(„_fc) j, h^" '^■'), j — 1,2, ... , dim(T-Of^), on the dual grid Di. Then 
the harmonic form h^^} becomes 

(4.36) *^(«.hr'^), with 

Note that, by construction, dh^^'^ = and d*h[^J = 0. 

4.4. Discussion. Although the reduction operator TZ and the reconstruction operator I are also 
introduced in [Sj, we have chosen to work in A^ = ITZK^ instead of = TZA^ . The reason 
is that metric concepts like an inner product, the wedge product and the Hodge operator are 
metric concepts which cannot be modeled in topology. Furthermore, these metric concepts are 
independent of the reduction operator and depend explicitly on the reconstruction operator I. 
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5. The Mimetic Spectral Element Method 

Now that a mimetic framework is formulated using differential geometry, algebraic topology 
and the relations between those, using the mimetic operators, we seek a projection that satisfies 
the properties of the mimetic framework. We restrict ourselves to piecewise polynomial recon- 
structions, as is common in most finite element/spectral element methods. 

First the polynomials with mimetic properties are described. Next the discretization is explained 
using different kinds of sample problems. 

5.1. Domain partitioning. Instead of solving in the infinite-dimensional space A'^(J7) we restrict 
our search to a proper subspace Af^i^^; Ck) C A'°(ri). Although there already exists several methods 
for the subspace A^(n;Cfe), (see for example |U [TH |^ ) , here we focus on a spectral element 
based method for curvilinear quadrilaterals. In spectral element methods a domain decomposition 
of fl is performed into M non-overlapping curvilinear quadrilateral closed sub-domains flm- 

M 

Q, = \^ fljn, Cl^l ~ dilrn H dill, m I, 

m—1 

where in each sub-domain a Gauss-Lobatto grid is constructed, see Figure [24j The collection of 
Gauss-Lobatto grids in all elements ftm constitutes the cell complex D. For each spectral element 
there exists a sub cell complex. Dm, i.e. the subcell complex in the spectral element flm- Note 
that Dm Di, m ^ I, is not an empty set in case they are neighboring elements, but contains all 
fc-cells, fc < n, of the common boundary, see Definition |30| 

Each sub-domain is obtained from the map : firef — >■ ^rm where the reference domain 
ilrcf is a unit cube flj-cf = [—1,1]", with n = dim(r2). All differential forms defined on Qm are pulled 
back onto this reference element using the following puUback operation, $*j : A''{il„i) — > A'^(i7i.ot). 
In three dimensions the reference element is given by 

(5.1) ^cf :={(C,^,C) I -1 <?,r/,C< 1}. 

The solution within each sub-domain is expanded with respect to a polynomial basis, corresponding 
to the chains in that element. 




X 



Figure 24. A 5 x 5 curvilinear multi-element Gauss-Lobatto grid. 

Remark 29. The map = ^^-1^ "is essentially the inverse of a chart from the computational 
manifold f2 to the unit cube in M". The map from the cells formed by the Gauss-Lobatto grid in 
[—1, 1]" to the cells in Q are the singular n-cubes. 

5.2. Reduction, mimetic basis-functions and projections. In the previous section four pro- 
jections - two direct projection and two coprojections - were introduced. In this section we derive 
mimetic basis-functions, that are cardinal basis-functions of arbitrary polynomial order, that are 
capable of reconstructing fc-cochains according to Definition [55j Because we consider only tensor 
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product based mimetic spectral elements for the reconstruction, it is sufficient to do the derivation 
and analysis in one dimension only. 

5.2.1. Projection tt^ using D. In spectral element methods the differential forms on Cl„i are approx- 
imated using piecewise polynomial expansions. The domain is given hy fl — Q^cf ■= ^ € [—1, 1]. 
We start with the projection operator tt/;, which is formed by the reconstruction of the reduction 
of a A;- form on the interior cell complex Di, see Section [3^ On il a cell complex D is defined ac- 
cording to Definition [30{ that consists of iV-f 1 nodes (0-cells), where — 1<Co<---<Cjv<1i 
and N line segments (1-cells), [^i-i, ^i], of which the nodes constitute the boundary, see Figure 20 
Consider a 0-form a'"^ — a(^) e A"(ri). Corresponding to this set of nodes (0-chain) there exists 
a projection, tt/i, using N^^ order Lagrange polynomials, li{£,), to approximate a 0-form, as 



(5.2) 



N 
i=0 



where ip is the isomorphism which associates a fc-c och ain to its coefficients in the expansion in 
terms of canonical basis A:-cochains, see Definition 
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The ith coefficient is referred to as tpi. 
Lagrange polynomials have the property that they interpolate nodal values and are therefore 
suitable to reconstruct the cochain a^^^ = 7?,a^''^(^) containing the set a, = a(^i) for i = 0, . . . , A^. 
So these polynomials can be used to reconstruct a 0-form from a 0-cochain. Lagrange polynomials 
are m fact 0-forms themselves, S A^(17rcf ; Co). They are constructed such that its value is 

one in the corresponding point and zero in all other grid points, 

(5.3) (Cp) '^'"^ 



if ij^p 



This satisfies (4.7), where in this case I = li'^^C)- It is straightforward to show that Lagrange 




Figure 25. The map from ^ e [^l, 1] hito a 1-dimensional sub-manifold in 



polynomials are invariant under a coordinate transformation. If li{x, y, z) is a Lagrange polynomial 
defined on a curvilinear 1-manifold embedded in a three-dimensional domain, then on that manifold 
there exist 0-cells T(o).p(a;, ?/, 2), associated to each node p, (xp,yp,Zp) , of the mesh of that manifold. 
In this case {x,y,z) = and {xp,yp,Zp) = ^(^p), see Figure 25 for a similar map in 2D. As a 
consequence ^*li{x,y, z) = li{S.), and so: 

Other useful properties of Lagrange polynomials are 



(5.4) 



TV 

E 

i=0 



TV 



z.(C) = i, 5]dZi(c) = o. 



i=0 
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Gerritsma [26 and Robidoux |70| independently derived the same projection for 1-forms, consisting 
of 1-cochains and l-form polynomials, that is called the edge polynomial, e\^\£,). 

application of the exterior derivative d to a^l^\S^), 



52 and 



Lemma 7. Following Definitions 
gives the l-form u'jl\S,) — dajj'^^) given by 



55 



N 



(5.5) 

with 1-cochain u^^\ where 
(5.6) 



daW(e) ^"^^ 



ni) 



a(0)(C) 



= al"^(C,)-«^"nC^-i)-fl.-a»-i, 
and with the edge polynomial defined as 

i-l N N i-1 

(0 = - E (a = E (0 - ^ E ) " ^ E (0 ■ 

fc=0 k=i k=i k=0 



(5.7) 



Proof. First take the exterior derivative of (|5.2 



d«r(o-E-'d^(°)(o. 

■i=0 

Now the l-form dal^-* is expanded in terms of the 0-cochain a*^"-' . Then a change of basis is applied 
to rewrite this expansion in terms of the 1-cochain Sa^'^'^ using the projection ttm, see Definition 61 
The trick is to subtract X^i^o "^^i^^lO: being equal to zero, and rewrite it such that we retrieve 
the coboundary operator and an edge polynomial, 

N 

= (0 = ^Mdar (e) ^ E(«^ - «^-)d4°^(0 

i=0 



N 

E 



E (aj-aj-i)+ E 



fc-i 



3 = k+l 
N 



d/r(c) 



Ed^^o E E d4°^(o E 



4=0 J=i+1 i=k+l 

k \ N 



j = k+l 



N 



E Ed^r^(^)h. + E Ed4°Ho 

j = l \i=0 / j=k+l 



N 



E-^(-E<'(o)*^E-^^i^^(^)- 



N 



fc=0 



□ 



The cochain corresponding to line segment (1-cell) T(i) ^ is given by Ui = ai — ai__i and so 
u'^) = (5a(°' is the discrete d eriva tive operator in ID. This operation is purely topolo gical; no 
metric is involved. It satisfies (4.8 1, since dXa^"^ = XSa'-^^\ Note that according to (2.11), we have 



dei(C) = dod^j^*^' (^) = 0. The polynomial l-form can be decomposed into a polynomial and the 
canonical basis for differential forms (see Proposition [3| , 



e['\0^e,m^, with e,(0--E^ 

k=0 ^ 
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Similar to (5.3), the edge functions are constructed such that when integrating e[^\^) over a hne 



segment it gives one for the corresponding element and zero for any other line segment, so 



(5.8) 



=(1) 



This also satisfies (4.7), where in this case I = e[^\^). The last property to verify is invariance 



1 if i — p 
if i^p 



under transformations. If ei{x,y,z) is an edge function defined on a curvilinear 1-manifold em- 
bedded in a three dimensional domain, then on that manifold there exist 1-cells r(i) p(x, y,z) 
associated to each edge, p, of the mesh of that manifold. In this case {x,y,z) — $(^), Figure 25 
and T(i)_p(a;, 2/, z) = ^{T^i-^iS,)). As a consequence ^*ei{x,y, z) — ei(^) and so: 



ei{x,y,z) 



(*-io$)(rf (0) 



This is what could be expected, since the generalized Stokes Theorem (2.15 1 is purely topological 



and does not depend on the particular coordinate system or polynomial representation. 

The following example shows the commutation property between projection and exterior deriv- 
ative. 



Example 20. In Figure 26 a graphical representation of the commutating diagram for the pro- 
jection and the exterior derivative is given. This figure illustrates Lemma [5[ 




Figure 26. An oscillatory 0-form (Top left) is projected onto a polynomial 0-form 
(Bottom left) using Lagrange functions. The exterior derivative of the oscillatory 
0-form (Top right) is projected onto polynomials using edge functions (Bottom 
right). This diagram commutes. 



Now that the mimetic basis-functions are defined it can be proven by example that the projec- 
tion, TTji, is in general not a Galerkin projection. 

Example 21. Let 11^ he a Galerkin projection, andYlf^a^^^ be expanded using edge basis-functions 
as Il/ja'"'^) — X^iLi^i^KO) then the coefficients di are determined by 

(n^a(i),e,(0) = {a<^'\e,iO). 
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5> 



GLL nodal interpolation 

I I I 



GLL edge interpolation 





Figure 27. Lagrange poly- 
nomials on Gauss-Lobatto- 
Legendre grid. 



Figure 28. Edge polynomials 
on Gauss-Lobatto-Legendre 
grid. 



In general, [ai ... un]"^ ^ TZa^^\ and therefore [ai 
a cochain projection. As an example, let a^^^ 



1 l]T 
4 ii ■ 



ajv] is not a cochain, and Hh is not 
x^dx, then for N = 1, [oi 0,2]"^ = i^lo 



Again, let N be the number of line segments in a spectral element. Then Lagrange basis 
interpolating N + 1 0-cells with fourth-order polynomials, corresponding to a Gauss-Lobatto grid, 
is shown in Figure [27] The edge basis interpolating A'^ line segments with third-order polynomials. 



corresponding to a Gauss-Lobatto grid, is shown in Figure 28 



5.2.2. Bounded linear projections. As mentioned in Definition |57t boundedness of the projection 
is a requirement, and is therefore shown for the projections introduced above. 

The mimetic framework uses Lagrange, li{^) G HA^{fl,-c{), and edge functions, ei(^) £ L^A^(rirof), 
for the reconstruction, I, where the latter is constructed using the former; i.e., from the finite di- 
mensional 0-form TT/ifl^^) = ^f^Qaili{^) G A° (fljcf ; Co), we define Tr/jfo^^' e A\{il.^cf',Ci), such 
that 

N 



da(°) 



Wjei(^). 



Because we consider tensor products to construct higher-dimensional interpolation, it is sufficient 
to show that the projection operator is bounded in one dimension. A similar approach was used 
in |14) . Due to the way the edge functions are constructed, there exists a commuting diagram 
property between projection and exterior derivative. 



HA" - 




> 








K - 




— > 0, 



which gives, for a(") G i7A"(r2,of), the 1-form 

(5.9) dW°) =7r,,da(°), in Ai(aef). 

Lagrange interpolation by itself does not guarantee a convergent approximation, [23j . but it re- 
quires a suitably chosen set of points, —1<^o<^i<-<£.n<^- Here, (extended)-Gauss and 
Gauss-Lobatto distributions are proposed, because of their superior convergence behaviour. The 
a priori error estimate for these kind of interpolants in the HA"-norm is given by, |16| 



(5.10) 



,(0) 



TThO 



(0) 



I if AO 



< C- 



p" 



/ min(p + 1, to). 
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Equation (5.101 also implies that the projection of zero- forms is bounded in the iJA°(f2ref); as is 
shown in the following proposition. 

Proposition 29. For a*-°-' € iJA°(f2i.cf) and the projection tt^ : HK^ — > A°, there exists the 
following two stability estimates in H -norm and H hP -semi-norm: 

(5.11) |k/,a(")||ffAo <C7||aW||^^Ao, 

(5.12) |W0)Uao <C|a(°)UAO. 

Proof. The iJA^'-norm of tt/jo'^''^ can be bounded from above using triangle inequality, 

For a(°) e iJA^, we use Poincare inequality, Lemmas) to bound the first term, 

||a("'lkAO <cp||da(0)|U.AO =cp|a(")UAO. 



Substituting this inequality and the interpolation estimate (5.101 into the triangle inequality gives 
the following result. 



(0) 



//AO- 



Then ( |5.11[ ) and \h.l2\ follow directly 

k/.a^°^|//Ao < Wa^''^\HAo < Cla^l^Ao < C\\a'^'^\\HAo. 



□ 



Now that we have a bounded linear projectiorj^of 0-forms in one dimension, we can also proof 
boundedness of the projection for 1-forms. 

Proposition 30. Let a'^^ G HA^ and b^^^ e i^A^, then the projection tt/i : L^A^ — > A^^ given by 
Lemma^is bounded 

(5.13) hhb^^^h^Ai <C||&(i)|U2Ai. 

Proof. Because L^A^ = 7?.(d; HA^), we can write b'^^^ = da'^^\ Then proof follows from the result 
of the previous proposition and the commutation between projection and derivative. Lemma [5] 

||7r„6(i)|U2Ai = |7r;,daWU.Ai - Id^.a^^^U^Ai = k^a^^^l/zAo 

< C|a(")|^A" = C|da(°)U.Ai = C||6<'^||l^ai. 

□ 

Propositions [29] and [30] show that the projection tt/i is a bounded linear projection, based on 
Lagrange functions and edge functions. Just like for 0-forms using Lagrange interpolation, we can 
also give an estimate for the interpolation error of 1-forms, interpolated using edge functions. 

Proposition 31. Let a(°) G HA° and 6(^' = da(°) G L^A\ the interpolation error 6^^' - 'Khb^^'> G 
L^A^ is given by 



(5.14) 

with I — min(p -|- 1, rri). 
Proof. 



- ^,.6(i)|U2Ai < C^|5(i)|^™-.Ai, 



and 



|aW-^,a(0)|^Ao<||a(")-W°'|kAo. 



Then (5.141 follows from (5.101, with the semi-norm rewritten as 



□ 



^also referred to as bounded cochain projection, [SJH]- 
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Corollary 11. Although Propositions \2S\ \30\ and \31\ were derived for the 0- and 1-forms on the 
reference domain fJicf . Due to the commuting property between the pullback and exterior derivative, 
Proposition^ these propositions also hold on any curvilinear domain Vl, where $ : iljcf — Q,. 

5.2.3. Projection tt/i using D. Given a cell complex D^, a corresponding dual grid Di was defined in 
Section [3. 3| For simplicity, consider a one-dimensional manifold, on which the two cell complexes 



are defined (e.g. Figure 20 1 . Let the Gauss-Lobatto grid defined above be the primal cell complex, 
Di, consisting of TV + 1 points and N line segments. The dual of D^, being Z)^, consists of N points 
and iV-l- 1 line segments, according to Definition |46j Having a Gauss-Lobatto primal cell complex, 
the points for the dual part could be the Gauss points, — 1 < < . . . < < 1, [IS]- We indicate 
the Gauss-Lagrange interpolant corresponding to Di by G A^(rii.cf ; Co), with Co{Di). The 
projection of a 0-form, a^°^ € A°{fl,.c{), using the 0-cells in D^, is a polynomial of degree N — 1, 
and is given by 



N 



(5.15) W°He) = E«»^"'f(^)- 

1=1 

The dual complex D is the union of Di and Dj,, where the two additional boundary points are 
given by = ^1 and £,n+i — 1- The 0-cells of the dual cell complex D are interpolated using 
extended-Gauss-Lagrange polynomials, ^^^(C) G ^hi^id'jCo), with Cq = Cq{D). The projection of 
a 0-form, a'-^' G A°(il), using the 0-cells in cell complex _D, is a polynomial of degree -I- 1, and 
is given by 

Af+l 

(5.16) = E 

■i=0 

The extended Gauss-edge functions, G K]^{Vl-ccf',Ci) with Ci G D, are found similarly to 
Lemma [t] The projection of a one-form b^^^ G A^(i7), using the 1-cells in D, is a polynomial of 
degree N , and is given by 

N+l i-1 

(5.17) ^,6W(0=E^'eT(0, with ~ef{(,)^-Y.dlf{0. 

i=l fc=0 

The extended-Gauss polynomials in the context of mimetic discretization were first preliminarily 
discussed in [101 El] . 

Remark 30. The boundary cells in Db can be seen as connectivity cells. Either they connect 
the computational domain to the outside world in terms of boundary conditions, or they connect 
adjacent spectral elements within the computational domain, thus providing C*^ -continuity for 0- 
forms. 

Remark 31. Although the boundary cells are part of the solution and thus should be solved for, 
the integral over the domain [—1, 1] of the boundary polynomials is zero, i.e. 

1 /.I 



1 j-i 

Therefore they do not contribute to the domain integrals overfl — [—1, 1]. There is a strong analogy 
with integration by parts. Let Qj^^ G A^^{Q]Cq) be represented using the primal cell complex and 
a"^^ G A^{Q]Co) be expressed using the dual cell complex. Then integration by parts gives 

d,r A 4"^ = - / A d4") + / A a("), ^ Kin; Co), 

n Jn Jan 

This shows that the derivative on the primal cell complex D is related to the derivative on the 
interior part of the dual complex, Di, supplemented with boundary part, Di,. 
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EGL nodal interpolation 




Figure 29. Lagrange poly- 
nomials on extended Gauss- 
Legendre grid. 



EGL edge interpolation 




Figure 30. Edge polynomials 
on extended Gauss-Legendre 
grid. 



Figure 29 shows the nodal basis functions on the dual grid, while Figure 30 shows the edge 
functions on the dual grid. Together with the basis functions depicted in Figures [27] and [28) 
these four sets of basis functions make up the entire set of reconstruction operators discussed in 
Section [H 



5.2.4. The coprojections and tt^. Now that the reconstruction polynomials on the primal and 
the dual complex have been introduced, the coprojections can be directly expressed in terms of 
these polynomials. For tt^, we take the Hodge-* of a fc-form a'^'^-' at the continuous level, which 
yields a (n — fc)-form. This form is projected by tt/j with respect to the dual grid. Then the 
Hodge-* operator is applied to the projected differential form and multiplied by (— This 
result is then expanded in terms of the /c-form basis functions on the primal grid. This is possible 



thanks to property (4.10) 



A similar route is followed by the coprojection tt^, where in this case the projection is with 
respect to the primal grid and the final result is expanded in terms of the fc-form basis functions 



on the dual complex. In Figure 31 the four projection operators are applied to a 0-form, a 



(0) 



sin(37ra:: -|- 0.08), and the resulting projected 0- forms arc plotted. 




Figure 31. Comparison between the four projections, tt, tt, it* and tt*, for a 
one-dimensional mesh of order N = 8. 
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5.3. Applications of discrete operators. In this section we will show some examples of the 
use of mimetic basis-functions and the action of the operators d, ★ and d* described in Section [4.2[ 
Now that Lagrange polynomials and edge polynomials are defined, we can construct interpolants 
for all fc-cochains in $7 C M". Because we consider quadrilateral elements only and employ tensor 
products to form the spectral element basis, the higher-dimensional basis functions are formed 
naturally by applying tensor products. For instance, in $7 C M'^, the surface element is the tensor 
product of two edge polynomials and one Lagrange polynomial, whereas the volume basis function 
is the tensor product of three edge polynomials. 

We start with the analog between the exterior derivative and the coboundary operator in 3D, 
as illustrated in the diagram below. 



n 



C°{D) 



-9> C\D) 



curl 



n 



div 



-9> C^D) 



We recognize the continuous and discrete gradient, curl and divergence operator from vector 
calculus. Again for clarity, we restrict ourselves to 51 = ^id, see (5.1). A similar De Rham 
complex can be set up for the dual complex using the nodal and edge functions given in the 
previous subsection. Note that in the finite dimensional setting there is no need to distinguish 



between and tt^, due to Proposition 24 



Example 22 (Gradient operator). Consider uj^^'' — dp^^\ where p^^'^ is expanded in standard 
coordinates (^, 77, C) as 



N N N 



(5.18) 



Ph 



(0) 



i=0 j=Q fe=0 

Apply the exterior derivative in the same way as in Lemma^ it gives 



(5.19) 



,(1) 



N N N 

i=i j=a k=o 

N N N 

i=0 j = l k=0 
N N N 

i=a j=a k=i 



where 
(5.20) 



and u[ 



ij,k ~ Phj,k Pi,j;k~l: 



can compactly be written as u'^^'^ = 6p^^\ or in matrix notation as 'ip{u^^')) — E^'^'^'>ilj{p^^'^). This 
relation is exact and invariant under transformations. Note that in terms of function spaces the 
vector proxies of p^^^ are in and that of are in _ff(curl). 
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Example 23 (Curl operator). Let u^j^^ be defined as in (5.19), then w'^^^ = dul^^-*. Apply the 
exterior derivative as in Lemma^ and consider the wedge product property (2.3), it gives 

N N N 

4'^ =EEE<.,/c^»(Oe,(ry)efe(C) 

2=0 j = l k=l 
N N N 

(5.21) +EEE<.,/ce^(O^.Wefe(C) 

i=l j=0 k=l 
N N N 

i=l 3 = 1 k=0 

where 

^ C C V _i_ V 

'^i,j,k ~ '^i,j,k ~ '^i,j-l,k ^ '^i,i,k ' "ij-./c-l' 
~ ""j.i,*: ~ '^i~l,j,k ~ "jj,fc + "ij-l,fe' 

can compactly be written as w^^^ = (Ju'-"'^), or m matrix notation as -(/'(w^^^) = E(^'"'^^'0(u("'^''). This 
relation is exact and invariant under transformations. Note that in terms of function spaces the 
vector proxies ofw'f^^ are in H{drv). Ifu\^^ is the gradient of p^l^\ then j j, — w^j j, — wjj ^ ^ 
and so uj^j^^ = 0. This is in accordance with 



'3.13) and (2.11). 



Example 24 (Divergence operator). Let the 2-form q'j^^ be defined as 

N N N 

i=0 j = l k=l 
N N N 

(5.23) +J2T.J2lLM^Mv)ekiC) 

1=1 j=0 k=l 
N N N 

+EEE4-,^e^(0e,(^)/fe(c), 

i=l j = l fe=0 

then u^'^' — dq'j^^ . Apply the exterior derivative as in Lemma and consider the wedge product 
property (2.3), it gives 

N N N 

(5.24) = E E E ^^^■.'=^^(^)^^ ('?)^fc(0. 

1=1 i=i fc=i 

where 

(5-25) Vij.k = qi^^k - 4-i^o,k + <llj,k - <ll]-i,k + 4.,j,k ~ 9i"j\fe-i' 

can compactly be written as v'^-* — 5ci^'^\ or in matrix notation as -ijjiy^^)) — E*^'^'^)'0(q*^^'). This 
relation is exact and invariant under transformations. It shows that dq'f^^ — dTr^g^^^ = dlTZq^^^ = 
XSTZq^'^^ = TTZdq^^^ = TT/jd^/^^^ — tt/jTj'"^) = W/f''. In case qj'^^ — w^j^^ = du^i^\ then it can be shown 
that ~ 0. Note that in terms of function .spaces the vector proxy of v'f^^^ is in L^ . 

Remark 32. The expansion of 0-forms a'j'^^ € A^(il;Co) has continuity over the sub-domain 
boundaries, dflm- For a^j^ G A^{Cl;Ci), the tangential component is continuous and for a^j^^ G 
A^j(i7; C2) the normal component hasC'^ continuity over the sub-domain boundaries. The expansion 
of volume- forms, a^(^'' G K\{D.;C^) are discontinuous over the sub-domain boundaries. 
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Example 25 (Trace). Without loss of generality, consider one- dimensional domain fl^^f ■ C € 
[—1,1]. We can apply the trace on a zero-form a^j^"^ G A^(r2i.cf ; Cq) and a one-form b^j^ £ 
A\{rt^ci',Ci), then we find: 



\i=0 
' N 



-ao, 



-On, 



N 



N 



' N 



'i=l 



tranH&i'^ =tran„ y^b,e,{^) = + ^fe,e,(+l). 

\i=l / i=\ 

The plus and minus signs in front indicate the corresponding orientation, see left in Figure So 
if we evaluate the trace on the left boundary point, where the orientation is negative, we obtain 
the positive value oq for the 0-form and + X^i^i f'^''" 1-form. Because we consider 

tensor products to consider higher- dimensional k-forms, the trace on k- forms in higher dimensional 
domains can be constructed straightforward using the above relations. 

Example 26 (Hodge-*). Now consider a ID domain fti-d : ^ — [—1, 1]. Let al^\^) = *'U'\p{£^). 
This example shows how to perform the Hodge-k. The action of the Hodge-k is followed by a 
projection ttm ■ A^{^rcf]Ci) ~> A^(r2i.of ; Co) to write the outcome in the preferred basis. Let 
'^h\0 expanded in terms of 1-cochains and edge- functions. 



,(1) 



N 



(0 = ^w,e,(0, with e Aiin;Ci). 



Then apply the Hodge-k to get a\^\s,) G A^(f2;Ci), as follows 



N 



N 



N 



Rewrite the basis using the action of 7Tj\/ 



,(0) 



N 

(0) 

'■h - ^Ma^ = 



N 



N 



In this case, what is usually called the Hodge-k matrix \30^. is given by (H'^'"'^)^^ = ei{£,j). The 
action of this Hodge-* is illustrated in Figure \3S\ 

Example 27 (Codifferential). Again consider a one- dim ension al domain f2].cf ■ C = [^li !]• ^6 
apply the codifferential on £ A^(Qj;c{;Ci) according to (2.28): 



N 

E 

J=0 



N 



fc=0 



The codifferential matrix D* is given by 
(5.26) (D*),_^. 



E 

fe=0 
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-1 -0.5 0.5 1 -1 -0.5 0.5 1 



f i 

Figure 32. Example of the action of the Hodge-* operator, fohowed by the 
projection ttm, t^m ° *■ Left the function and corresponding cochain for u^^^^ g 
A^(ri; Ci). Right the function and corresponding cochain for af^'^ G A"(J7; Co). 



The coefficients in this codifferential matrix are independent of the location of the dual grid, as 
would be expected. Even when the codifferential matrix is constructed using two Hodge matrices and 
an incidence matrix, identically the same codifferential matrix is found. The same is true when 



retrieving the codifferential matrix from the formal adjoint formulation {2.26) when the 1-form 



kdif formhul is zero at the boundary. Note that (5.26) requires that li{^) is twice differ entiable 



while for the adjoint formulation it is sufficient that li{C) is only one time differ entiable. 

Example 28 (Hodge decomposition). In this example we show the result of the discrete Hodge 
decomposition applied to a velocity field, Vjp, obtained from a potential flow problem on an annulus. 
The Hodge decomposition of the velocity field is given in terms of the gradient of the potential, 
and a harmonic function, h^j^\ 

For the annulus, consider the domain {r,9) € x [0, 27r] and the cell complex D which covers 

this domain as shown in Figure The topology of the cell complex is equal to the topology of the 
cell complex shown in Figure U 




Figure 33. The cell complex which covers the {r,9) e [l,R] x [0,2n] and the 
1-cochain values associated with each 1-cell. 
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(5.27) 



The velocity 1-form on the annulus is given by 

1 



.,(1) 



1 



dr 



1 



r 

27rr 



de 



Application of the reduction map, TZ, gives the 1-cochain values associated with each 1-chain in 
the cell complex. Its values are given in Figure If we apply the cohoundary to this 1-cochain 
we get the zero 2-chain. Because the topology is the same as in Figure the harmonic 1-chain 
is still given by 

h(i) = (1,-1,0,0,1,1,-1,1,-1,0,0,-1)^ , 

and therefore 

r 



a 



a — 

4 



The harmonic cochain h^^J_was determined in Example IS The values of the 1-cochain 
^h^"'^^ are given in Figure 33 The projection of the velocity field becomes, 



s(o) 



TThV 



(1) 



I 5cj) 



(0) 



.h(i) 



Figure \3Ji\ shows the reconstruction of the velocity field components and the total velocity field, for 
different vortical strengths and different polynomial order. 



Non-lifting flow 

2 p = 4 





Figure 34. Decomposition of velocity field. It shows the curl-free components 
(left), the harmonic components (top) and the total velocity field, for different 
vortical strengths and different polynomial order. 
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5.4. Discussion. Closely related to the presented mimetic basis-functions are (higher-order) Whit- 
ney forms, [511211122]. They share that there exists points, edges, faces and volumes (0,1,2,3-cells) 
both on the element boundary as well as in the element interior. These fc-cells have physical 
relevance [66]. The main difference exists in the way they are constructed. 

On the other hand in mixed finite elements, Raviart-Thomas [67' (iJ(div) conforming) and Nedelec 
[?3j (-ff (curl) conforming) elements are frequently used basis- functions. Main difference with our 
mimetic basis-functions is that these basis- functions are moment based polynomials, where each k- 
cell may have a higher-order moment degree, but the number of /c-cells remains unchanged when 
increasing the polynomial order. A consequence is that Raviart-Thomas and Nedelec elements 
allow only affine mappings, whereas the present mimetic basis-functions allow for curvilinear ele- 
ments. 

Where the mimetic basis-functions are constructed especially for quadrilateral and hexahedral 
elements, the Whitney and mixed finite element basis-functions also work on triangles and tetra- 
hedrals. In case of the lowest order quadrilaterals or hexahedrals, all three types of basis-functions 
are the same. 

As shown in Section|3] we need two cell-complexes to represent inner- and outer oriented variables. 
The use of such overlapping grids is well-known in staggered finite volume methods. A spectral 
staggered grid approach was presented by Kopriva and Kolias, [50], and Kopriva, [38t i39j . 
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6. Coda 

Topological notions in PDE's. In the mimetic framework we presented, building bloclcs 
are introduced, which mimic in a finite dimensional setting differential operators, which appear 
in PDE's as used in physics. By mimicking we mean that operators act similarly in a continuous 
space. A*-', to the operations in a finite dimensional space, A^, or the discrete space, . This is 
refiected by the many commuting diagrams presented throughout this paper. The representation 
of a finite dimensional solution is directly related to the underlying mesh. We made a clear 
distinction between topological relations that are metric free and exact, and metric relations that 
contain the constitutive relations and approximations. We extensively described the topology 
of the mesh, consisting of fc-cells that have either inner or outer orientation, the basis of a cell 
complex. Topology and orientation are intrinsically related to differentiation as is indicated by 
the boundary and coboundary operators and their corresponding complexes. Key ingredient for 
exact differentiation is the generalized Stokes' Theorem. This idea is not new and is commonly 
referred to as finite volume methods, where the discrete unknowns, i.e. the cochains, are integral 
quantities associated to geometric objects, obtained by the reduction map. In this way the action 
of differential operators is independent of the reconstruction map. Illustrations of this fact were 
given in Examples [22l [23l and [^1 

Metric concepts in PDE's. The existence of constitutive relations necessitates us to relate 
physical quantities associated to geometric objects of different dimension and different orientation. 
The one-to-one relation between physical quantities connected by constitutive relations motivates 
the existence of a dual grid, as is known from (staggered) finite volume methods. In dual grid 
methods, these metric relations are usually treated by a finite dimensional representation of the 
Hodge-* operator. This finite dimensional Hodge-* is not unique, but depends on the choice of 
reconstruction operator, in our case by arbitrary-order mimetic spectral element interpolation. In 
contrast, finite element methods usually consider i^-inner products as metric operator. The two 
are in fact intrinsically related using the wedge product according to Definition [20| While in finite 
volume reconstruction is usually seen as interpolation, in finite elements the reconstruction func- 
tions are better known as basis-functions, ignoring the fact that these 'functions' are differential 
forms which make the connection with the associated geometric objects. 

Single grid methods. In terms of the framework presented in this paper, in the finite element 
method the weighting function 'lives' on the same cell complex as the complex in which the 
equation is defined. The Hodge operator takes the weighting functions to the dual complex and 
then the wedge product is evaluated. So, although in finite element methods, one generally does 
not construct a dual cell complex, the dual cell complex is implicitly incorporated through the 
inner product. When we recognize that the weighting functions refer to the dual cell complex, it 
is also possible to make a more physical interpretation of integration by parts in finite element 
methods. Integration by parts transfers an equation from the dual complex to the primal complex 
(and the other way around). 

Dual boundaries. Orientation not only introduces primal and dual cell complexes, it intro- 
duces primal and dual grids for both the domain and the boundary dVL. This distinction is 
directly related to integration by parts. Especially for the dual complex, the boundary dual com- 
plex is often not recognized in staggered finite volumes methods, but instead additional degrees 
of freedom are generated ad hoc which are called 'ghost points'. 

Hodge decomposition. This framework presented more than just a nice analogy between 
finite volume and finite element methods. It has more structure. Most important is the Hodge 
decomposition and its finite dimensional and discrete counterparts. The Hodge decomposition 
decomposes differential form spaces into separate function spaces related to the exterior derivative. 
Mimicking the Hodge decomposition is essential for numerical stability. A key ingredient here is 
the commutation between projection operator and the exterior derivative. 

Extension to curvilinear grids. Finally it has been shown that all relations hold for both 
Cartesian, as well as curvilinear grids, thanks to the preservation of the commutation relations of 
the pullback and cochain map with the various operators discussed in this paper. The topological 
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relations are insensitive to changes in the topology as long as the grid connectivity remains unal- 
tered, whereas the metric-dependent part (inner-product and wedge) do change when one maps 
the standard element to a curvilinear domain. 
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